root / src / Path.f90 @ 8
History  View  Annotate  Download (83.7 kB)
1 
! This programs generates and optimizes a reaction path 

2 
! between at least two endpoints. 
3 
! 
4 
! It uses NAMELIST for input, see below ! 
5 
! 
6 
! v3.0 
7 
! First version entirely in Fortran. 
8 
! It uses routines from Cart (path creation and parameterisation) 
9 
! and from IRC06, especially routines for point optimization. 
10 
! 
11 
! TO DO: 
12 
! 1) Possibility to read and/or calculate some hessian structures 
13 
! 2) Use of internal coordinate to estimate the hessian for stable 
14 
! species in a much better way... 
15 
! 
16 
!!!!!!!!!!!!!!! 
17 
! v3.1 
18 
! The Hessian update includes neighbour points. 
19 
! 
20 
! v3.2 
21 
! We add the option Hinv that uses the BFGS update of the Hessian inverse 
22 
! 
23 
! v3.3 
24 
! We add the full option for ZMAT coordinates, where all is done using 
25 
! internal coordinates. 
26 
! v3.31 
27 
! The step size is controlled using the step norm and not the maximum step 
28 
! component. 
29 
! v3.5 
30 
! Big modifications of Calc_zmat_const_frag in order to be able 
31 
! to treat VASP geometries. 
32 
! For now, only XYZ geometries are read and written. VASP interface 
33 
! has to be written. 
34 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
35 
! v3.6 
36 
! Some minor modif in Calc_zmat_const_frag to have nicer programming 
37 
! and a new mode added: Coord=MIXED where part of the system 
38 
! is extrapolated in cartesian coordinates and the rest of it in zmat. 
39 
! it might be Useful for VASP, or when lots of atoms are frozen 
40 
! by default, when Coord=MIXED, all frozen atoms are described in cartesian 
41 
! coordinates. 
42 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
43 
! v3.61 
44 
! We move further and further to VASP program. New variables added to the 'path' namelist: 
45 
! input: what is the program used for the geometries ? 
46 
! Possible values: xyz (or cart) for Xmol format; VASP. Xyz should be used for Gaussian 
47 
! prog: what is the program used for the calculations ? 
48 
! Possible values for now: gaussian, vasp. 
49 
! For now (02/2007), when prog='vasp' is used, only the initial path is created. 
50 
! That means that we have only added input/output subroutines :) 
51 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
52 
! v3.7 
53 
! We try to improve the optimization steps. For that, we first modify many routines 
54 
! so that the new step is calculated in the 'free' degrees of freedom, as it is done 
55 
! in IRC07/ADF for example. That will allow us to impose constraints in an easier way. 
56 
! 
57 
! v3.701 
58 
! Add a new option for the initial extrapolation of the Hessian 
59 
! 
60 
!v3.71 
61 
! The optimization step has been improved... but the vfree part has not yet been included. 
62 
! This has still to be done :) 
63 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
64 
! v3.8 
65 
! I have included part of the vfree procedure: for now freemv returns the Id 
66 
! matrix, that is used to project out the tangent vector. 
67 
! This improves the optimization a lot, but I will have to implement 
68 
! a real freemw subroutine to have constraints there ! 
69 
! v3.81 Technical version 
70 
! I have added a new program 'Ext'. When this one is used, Path creates a Ext.in file, 
71 
! calls Ext.exe file to get a Ext.out file that contains the Energy and the gradient. 
72 
! Format: 
73 
! Ext.in : Xmol format. The second line (comment in Xmol) contains the COORD name. 
74 
! Ext.out: Energy on firt line (1X,D25.15), then gradient in CARTESIAN coordinates 
75 
! (3(1X,D25.15)) format. Natoms lines. 
76 
! TO DO: Gradient could be given in COORD format, ie cartesian for COORD=CART, XYZ or HYBRID 
77 
! ZMAT for CORD=ZMAT (in that case, 6 zeros are there at the begining !) 
78 
! For MIXED, it is equivalent to CART :) 
79 
! v3.811 
80 
! I have added a prog="TEST" option, where egradient calls egrad_test subroutine. 
81 
! It then calculates the energy and gradient as it wants and it returns the cartesian 
82 
! gradient. The purpose is to have an analytical approximation of the PES of the system 
83 
! under study to gain some time in testing the program. This is not documented :) 
84 
! 
85 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
86 
! v3.9 
87 
! In order to be fully interfaced with VASP, I have to change the architecture 
88 
! of the program. In particular, with VASP, one can use a Para8+NEB routine 
89 
! of VASP to calculate energy+forces for all points of the path. 
90 
! v3.91 
91 
! One small modification: in the cart list, one can freeze a whole sequence by 
92 
! entering a negative number for the last atom of the sequence 
93 
! for example : cart=1 160 165 170 will define atoms from 1 to 160 and from 
94 
! 165 to 170 (included) as cartesian atoms. 
95 
! Of course, Same applies to Frozen 
96 
! v3.92 / v3.93 
97 
! Slight modifications. One noticeable: when using Vasp, the program analyses the 
98 
! displacements to suggest which atoms could be described in CART. 
99 
! v3.94 
100 
! VaspMode=Para now implemented. One needs to be careful in defining 'ProgExe'. 
101 
! In particular, one must put the mpirun command in ProgExe !!! There is NO test. 
102 
! v3.95 
103 
! When using internal coordinates (zmat, hybrid or mixes), Path calculates the number 
104 
! of fragments for each geometry. The reference one (ie the one used to compute the internal 
105 
! coordinates) is now the one with the least fragments. 
106 
! user can also choose one geometry by specifying IGeomRef in the Path namelist. 
107 
! v3.96 
108 
! I have added a test for the geometry step: it does not allow valence angle to 
109 
! be negative or bigger than 180. I hope that this might help geometry optimization 
110 
! when dealing with nearly linear molecules. 
111 
! v3.961 
112 
! The options cart and frozen are now LOGICAL: 
113 
! Fcart=T indicates that a &cartlist namelist should be read ! 
114 
! Ffrozen=T indicates that a &frozenlist namelist should be read ! 
115 
! TO DO: Autocart=T indicates that the atoms described in cartesian coordiantes should 
116 
! be found automatically 
117 
! v3.97 
118 
! Baker has been implemented by P. Dayal: we are interpolating Xint for now. 
119 
! TO do: interpolate Umat ? 
120 
!  
121 
! v3.971 
122 
! AutoCart=T is being implemented. Goal: having a 'black box' version for Vasp users. 
123 
! Vmd: True if user wants to watch the path using VMD. Only used for VASP for now. 
124 
! 
125 
! v3.98 
126 
! Calc_zmat has been replaced by Calc_zmat_frag which is more robust. 
127 
! Still missing: linear molecules and some tests. 
128 
! 
129 
! v3.981 
130 
! * Added a non document flag in .valid file that switch on the numerical 
131 
! calculation of the hessian instead of the update. 
132 
! * Added HesUpd variable in the namelist. Default value is MS for path 
133 
! optimization because it does not imposes a positive hessian, and BFGS 
134 
! for geometry optimization 
135 
! v 4.0 
136 
! * This version has been made public within the CARTE project. 
137 
! * Added: 
138 
!  Step_method to choose the optimization method: RFO or Diis 
139 
!  DynMaxStep: if TRUE, the maximum size of a step is updated at each step. 
140 
! If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. 
141 
! It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0] 
142 
! 
143 
!!!!!!!!!!!!!!!!!!!!!! 
144 
! v4.1 
145 
! * Para mode has been partly implemented for Gaussian. 
146 
! VaspMode is thus now 'RunMode' and can be Serial or Para 
147 
! Only Gaussian and VASP are supported for Para mode now. 
148 
! v4.12 
149 
! Some bugs in unconstrained zmat construction are corrected. 
150 
! v4.13 
151 
! Prakash has implemented the GEDIIS optimization procedure. 
152 
! v4.14 
153 
! There is a small problem for some inputs in VASP, where the 
154 
! last geometry is not the one given by the user. 
155 
! It seems to come from the fact that some people try to freeze 
156 
! some unfrozen atoms 
157 
! So some tests have been added to check that frozen atoms do not move. 
158 
!!! 
159 
! v4.15 
160 
! There were some bugs when reading and expanding cartlist and frozenlist 
161 
! I have changed the way frozen atoms are compared by adding partial alignment. 
162 
! v4.16 
163 
! Some bugs were corrected. 
164 
! v4.17 
165 
! Added MOPAC as a program. 
166 
! v4.175 
167 
! Added Three more parameters for defining the input and output names for 
168 
! the calculations. 
169 
! CalcName: Prefix for the files used for the energy and gradient calculations 
170 
! ISuffix: Suffix for the input file 
171 
! OSuffix: suffix for the output file. 
172 
! This maters for Gaussian for example. 
173 
! 
174 
! v4.177 (P) 2010 Nov 22 
175 
!  We add TurboMole as a new external code. 
176 
!  Subtle change for Input: the default is XYZ whatever the code, except 
177 
! for VASP. 
178 
! 
179 
! v4.178 (P) 2010 Nov 28 
180 
!  We add the possibility to call CHAMFRE reactive force field using 
181 
! prog=TEST2 or CHAMFRE 
182 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
183 
! OpenPath v1.4 2011 MayJune 
184 
!  Also added two flags: CalcEReac and CalcEProd. If .TRUE. Path will 
185 
! compute the energy of the reactant 
186 
! (and/or Product) at the begining of the calculations. 
187 
! This is useful if you start from non optimized points. 
188 
! 
189 
!  Also added more methods for Hessian updates. The following methods are 
190 
! now availables: 
191 
! For inverse Hessian: BFGS, DFP, MS. 
192 
! For Hessian: BFGS, DFP, MS, Bofill 
193 
! 
194 
! Small change: HesUpd is replaced by Hupdate 
195 
! 
196 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
197 
! TO DO 
198 
!  Trying to improve the output of Path. In particular, we now create 
199 
! a .datl file containing the image index, curvilinear coordinate and energy 
200 
! in the spirit of the AnaPath tool. 
201 
! AnaPath is still usefull as it can extract geometrical parameters on top 
202 
! of the energy 
203 
!  When using Hessian update, it is more efficient to use Cholesky 
204 
! decomposition, and to update this decomposition rather than the Hessian. 
205 
! See Nocedal&Wright p141. 
206 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
207 
! 
208 
! Project Name has changed to Opt'n Path 
209 
! OpenPaht 2011 MayJune > Opt'n Path v1.4 
210 
! 
211 
! 2013 Jan 18 Opt'n Path v1.47 (going to 1.5 slowly) 
212 
! We add Siesta (and start thinking of add CP2K) 
213 
! as a new "energy engine" code 
214 
! TO DO to go to 1.5: finish cleaning... 
215 
! 
216 
! 2013 Feb Opt'n Path v1.48 
217 
! We add some keyword for more output for Geometry Optimizations 
218 
! GeomFile: name of the file to print out the geometries and their energy 
219 
! as Xyz file. (only format for now) 
220 
! 
221 
! 2013 Feb opt'n Path v1.49 
222 
! We add the possibility to analyze the geometries in terms 
223 
! of Bonds, Angles, dihedrals 
224 
! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path) 
225 
! This is basicaly the merging of AnaPath into Opt'n Path... at last ! 
226 
! This is done by adding a Flag into the Path namelist: 
227 
! AnaGeom: logical, if true Path looks for the AnaList namelist 
228 
! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path 
229 
! will read it and print the values of the variable in a Path.dat file 
230 
! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. 
231  
232 
PROGRAM PathOpt 
233  
234 
use VarTypes 
235 
use Path_module 
236 
use Io_module 
237  
238 
IMPLICIT NONE 
239  
240 
CHARACTER(132) :: FileIn, FileOut, Line 
241  
242 
INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF 
243 
INTEGER(KINT) :: IOpt 
244 
INTEGER(KINT) :: Idx 
245 
INTEGER(KINT) :: N3at, NFree, Nfr 
246  
247 
INTEGER(KINT) :: List(2*MaxFroz) 
248  
249 
REAL(KREAL) :: E,FactStep !E=Energy 
250 
REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad 
251 
REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E 
252 
REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord 
253 
REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord) 
254 
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord 
255 
REAL(KREAL), ALLOCATABLE :: HessTmp(:,:) 
256 
REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim) 
257 
REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat6)) 
258 
LOGICAL, SAVE :: allocation_flag 
259  
260 
! Flag to see if frozen atoms are frozen (see below) 
261 
LOGICAL :: FChkFrozen 
262  
263 
! Energies for all points along the path at the previous iteration 
264 
REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF 
265 
! Maximum step allowed for this geometry 
266 
REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF 
267  
268 
! these are used to read temporary coordinates 
269 
LOGICAL :: FFrozen,FCart 
270  
271 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
272 
! 
273 
! To analyse the frozen atoms 
274 
! 
275 
REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat 
276 
REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat 
277 
REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist 
278 
LOGICAL, ALLOCATABLE :: ListAtFroz(:) 
279 
INTEGER(KINT) :: Iat1, Iat2, Iat3 
280  
281 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
282  
283 
LOGICAL :: Debug, Fini,Flag_Opt_Geom 
284  
285 
! INTEGER(KINT), EXTERNAL :: Iargc 
286 
INTEGER(KINT), EXTERNAL :: ConvertNumAt 
287 
INTEGER(KINT) :: IoS 
288  
289 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
290 
! 
291 
! For Debugging purposes 
292 
! 
293 
LOGICAL :: FCalcHess 
294 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
295  
296  
297 
INTERFACE 
298 
function valid(string) result (isValid) 
299 
CHARACTER(*), intent(in) :: string 
300 
logical :: isValid 
301 
END function VALID 
302  
303 
SUBROUTINE Read_Geom(Input) 
304 
CHARACTER(32), INTENT(IN) :: input 
305 
END SUBROUTINE READ_GEOM 
306  
307  
308 
subroutine hupdate_all (n, dx, dgrad, HessOld) 
309  
310 
IMPLICIT NONE 
311  
312 
INTEGER, PARAMETER :: KINT=KIND(1) 
313 
INTEGER, PARAMETER :: KREAL=KIND(1.0D0) 
314  
315 
INTEGER(KINT), INTENT(IN) :: n 
316 
real(KREAL) :: dx(n), dgrad(n), HessOld(n*n) 
317 
END subroutine hupdate_all 
318  
319 
SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess) 
320  
321 
IMPLICIT NONE 
322 

323 
INTEGER, PARAMETER :: KINT=KIND(1) 
324 
INTEGER, PARAMETER :: KREAL=KIND(1.D0) 
325 

326 
INTEGER(KINT), INTENT(IN) :: NFREE 
327 
REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree) 
328 
REAL(KREAL), INTENT(IN) :: ds(Nfree) 
329 
REAL(KREAL), INTENT(IN) :: dGrad(Nfree) 
330 

331 
END subroutine Hinvup_BFGS_new 
332  
333  
334 
FUNCTION InString(Line,String,Case,Clean,Back) Result(Pos) 
335  
336 
Use VarTypes 
337  
338 
implicit none 
339 
! Input 
340 
CHARACTER(*), INTENT(IN) :: Line 
341 
CHARACTER(*), INTENT(IN) :: String 
342 
LOGICAL, OPTIONAL, INTENT(IN) :: Case 
343 
LOGICAL, OPTIONAL, INTENT(IN) :: Back 
344 
CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean 
345  
346 
! Output 
347 
! the position of String in Line (the first one) unless Back is present 
348 
INTEGER(KINT) :: Pos 
349 
END FUNCTION InString 
350  
351 
SUBROUTINE die(routine, msg, file, line, unit) 
352  
353 
Use VarTypes 
354 
Use io_module 
355  
356 
implicit none 
357 
character(len=*), intent(in) :: routine, msg 
358 
character(len=*), intent(in), optional :: file 
359 
integer(KINT), intent(in), optional :: line, unit 
360  
361 
END SUBROUTINE die 
362  
363 
SUBROUTINE Warning(routine, msg, file, line, unit) 
364  
365 
Use VarTypes 
366 
Use io_module 
367  
368 
implicit none 
369  
370 
character(len=*), intent(in) :: routine, msg 
371 
character(len=*), intent(in), optional :: file 
372 
integer(KINT), intent(in), optional :: line, unit 
373  
374 
END SUBROUTINE Warning 
375  
376  
377 
END INTERFACE 
378  
379 
Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, & 
380 
MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, & 
381 
SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, & 
382 
input,prog, ProgExe,RunMode, AtTypes,poscar,PathName, & 
383 
OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, & 
384 
FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, & 
385 
NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, & 
386 
NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile 
387  
388 
Namelist/cartlist/list 
389 
Namelist/frozenlist/list 
390 
Namelist/analist/nb 
391  
392  
393 
Flag_Opt_Geom = .FALSE. ! Initialized. 
394  
395 
!!!!!!!!!!!!!!!!!!!!!!!!! 
396 
! 
397 
! We print the Version info in file 
398 
! 
399 
WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 20072013" 
400  
401 
! We read the files names 
402 
! SELECT CASE (iargc()) 
403 
SELECT CASE (command_argument_count()) 
404 
CASE (2) 
405 
call getarg(1,FileIn) 
406 
OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') 
407 
call getarg(2,FileOut) 
408 
OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN') 
409 
CASE (1) 
410 
call getarg(1,FileIn) 
411 
IF ((FileIn=="help").OR.(FileIn=="help")) THEN 
412 
WRITE(*,*) "Path minihelp" 
413 
WRITE(*,*) "" 
414 
WRITE(*,*) "" 
415 
WRITE(*,*) "Use: Path Input_file Output_file" 
416 
WRITE(*,*) "Input_file starts with a Namelist called path" 
417 
WRITE(*,*) "" 
418 
WRITE(*,*) "Compulsory variables are:" 
419 
WRITE(*,*) "NGeomi: Number of geometries defining the Initial path" 
420 
WRITE(*,*) "NGeomf: Number of geometries defining the Final path" 
421 
WRITE(*,*) "Nat : Number of atoms" 
422 
WRITE(*,*) "" 
423 
WRITE(*,*) "" 
424 
WRITE(*,*) "Other options" 
425 
WRITE(*,*) "Input: string that indicates the type of the input geometries." 
426 
WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp" 
427 
WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations." 
428 
WRITE(*,*) " Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext" 
429 
WRITE(*,*) " * In case of a Gaussian or Mopac calculations, input must be set to Cart." 
430 
WRITE(*,*) " One example of a gaussian/mopac input should be added at the end of the" 
431 
WRITE(*,*) " input file. See example file Test_G03.path or Test_Mopac.path" 
432 
WRITE(*,*) " * In the case of a VASP calculation, if input is set to Cart, then" 
433 
WRITE(*,*) " the preamble of a VASP calculation should be added at the end of" 
434 
WRITE(*,*) " the input file. See example file Test_VASP_cart.path" 
435 
WRITE(*,*) " In the case of a VASP calculation, one should also give value of the " 
436 
WRITE(*,*) " Runmode variable" 
437 
WRITE(*,*) " * In the case of a SIESTA calculation, an example of a Siesta input file" 
438 
WRITE(*,*) " should be added at the end of the input file. See Test_Siesta.path" 
439 
WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy" 
440 
WRITE(*,*) " and gradient of the whole path or not. If one wants to use VASP," 
441 
WRITE(*,*) " Runmode must be set to PARA." 
442 
WRITE(*,*) "PathName: Prefix used to save the path. Default is Path" 
443 
WRITE(*,*) "Poscar: string that will be used as the prefix for the different " 
444 
WRITE(*,*) " POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.," 
445 
WRITE(*,*) " not used for internal calculations." 
446 
WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations" 
447 
WRITE(*,*) " ISuffix: Suffix for the input file." 
448 
WRITE(*,*) " OSuffix: suffix for the output file." 
449 
WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates." 
450 
WRITE(*,*) " Valid only for Coord=Zmat, Hybrid or Mixed" 
451 
WRITE(*,*) "Fact: REAL used to define if two atoms are linked." 
452 
WRITE(*,*) " If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked." 
453 
WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info." 
454 
WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:" 
455 
WRITE(*,*) "  CART (or Xyz): works in cartesian" 
456 
WRITE(*,*) "  Zmat: works in internal coordinates (Zmat)" 
457 
WRITE(*,*) "  Mixed: frozen atoms, as well as atoms defined by the " 
458 
WRITE(*,*) " 'cart' array(see below) are describe in CARTESIAN, whereas the" 
459 
WRITE(*,*) " others are described in Zmat" 
460 
WRITE(*,*) "  Baker: use of Baker coordinates, also called delocalized internal coordinates" 
461 
WRITE(*,*) "  Hybrid: geometries are described in zmat, but the gradient are used in cartesian" 
462 
WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
463 
WRITE(*,*) "  RFO: Rational function optimization" 
464 
WRITE(*,*) "  GDIIS: Geometry optimization using direct inversion in the iterative supspace" 
465 
! WRITE(*,*) " HesUpd: Deprecated. Use Hupdate." 
466 
WRITE(*,*) " HUpdate: method to update the hessian. By default, it is MurtaghSargent" 
467 
WRITE(*,*) " Except for geometry optimization where it is BFGS." 
468 
WRITE(*,*) " GeomFile: Name of the file to print the geometries and their energy" 
469 
WRITE(*,*) " during a geometry optimization. If this variable is not given" 
470 
WRITE(*,*) " then nothing is printed." 
471 
WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization" 
472 
WRITE(*,*) "Smax: Maximum length of a step during path optimization" 
473 
WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary" 
474 
WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken" 
475 
WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path" 
476 
WRITE(*,*) " that is kept. FTan=1. corresponds to the full displacement, " 
477 
WRITE(*,*) " whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0." 
478 
WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization." 
479 
WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path." 
480 
WRITE(*,*) " This option indicates the first step where spline interpolation is used." 
481 
WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, " 
482 
WRITE(*,*) " it might happen that an atom moves out of the unit cell." 
483 
WRITE(*,*) " Path detects this by comparing the displacement to boxtol:" 
484 
WRITE(*,*) " if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5" 
485 
WRITE(*,*) "" 
486 
WRITE(*,*) "Arrays:" 
487 
WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements." 
488 
WRITE(*,*) " You can modify it using, rcov(6)=0.8." 
489 
WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements." 
490 
WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations." 
491 
WRITe(*,*) "If not given, Path will read the POTCAR file." 
492 
WRITE(*,*) "" 
493 
WRITE(*,*) "Flags:" 
494 
WRITE(*,*) "MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE." 
495 
WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time." 
496 
WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products." 
497 
WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants." 
498 
WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para" 
499 
WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para" 
500 
WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
501 
WRITE(*,*) "Hinv: if True, then Hessian inversed is used." 
502 
WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations." 
503 
WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path." 
504 
WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
505 
WRITE(*,*) " If True, a &frozenlist namelist containing the list of frozen atoms must be given." 
506 
WRITE(*,*) " If VASP is used, and frozen is not given" 
507 
WRITE(*,*) " here, the program will use the F flags of the input geometry" 
508 
WRITE(*,*) "FCart: True if one wants to describe some atoms using cartesian coordinates. " 
509 
WRITE(*,*) " *** Only used in 'mixed' calculations. ***" 
510 
WRITE(*,*) " If True, a &cartlist namelist containing the list of cart atoms must be given." 
511 
WRITE(*,*) " By default, only frozen atoms are described in cartesian coordinates." 
512 
WRITE(*,*) "" 
513 
WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis." 
514 
WRITE(*,*) " If True, Opt'n Path will look for the AnaList namelist after the Path Namelist." 
515 
WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then" 
516 
WRITE(*,*) " Opt'n Path will read it and print the values of the variable in a .dat file" 
517 
WRITE(*,*) " If AnaGeom is T but Analist is not given, then Path.dat just contains the energy." 
518 
WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description" 
519 
WRITE(*,*) "of the variables among:" 
520 
WRITE(*,*) "b(ond) At1 At2" 
521 
WRITE(*,*) "a(ngle) At1 At2 At3" 
522 
WRITE(*,*) "d(ihedral) At1 At2 At3 At4" 
523 
WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added" 
524 
WRITE(*,*) "at the end of the real atoms of the system" 
525 
WRITE(*,*) "" 
526 
WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry." 
527 
WRITE(*,*) " If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. " 
528 
WRITE(*,*) " It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]" 
529 
WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms." 
530 
WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now" 
531 
WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files." 
532 
WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP." 
533 
WRITE(*,*) "" 
534 
STOP 
535 
END IF 
536  
537 
OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') 
538 
IOOUT=6 
539 
CASE (0) 
540 
IOIN=5 
541 
IOOUT=6 
542 
END SELECT 
543  
544  
545 
! We initiliaze variables 
546 
Pi=dacos(1.0d0) 
547 
Nat=1 
548 
Fact=1.1 
549 
IGeomRef=1 
550 
NGeomI=1 
551 
NGeomF=8 
552 
IReparam=5 
553 
IReparamT=1 
554 
ISpline=5 
555 
debug=.False. 
556 
MW=.TRUE. 
557 
bohr=.False. 
558 
Coord='ZMAT' 
559 
Step_method='RFO' 
560 
DebugFile='Path.valid' 
561 
PathName="Path" 
562 
MaxCyc=50 
563 
SMax=0.1D0 
564 
SThresh=1. 
565 
GThresh=1. 
566 
FTan=0. 
567 
Hinv=.TRUE. 
568 
IniHUp=.FALSE. 
569 
HupNeighbour=.FALSE. 
570 
HesUpd="EMPTY" 
571 
HUpdate="EMPTY" 
572 
FFrozen=.FALSE. 
573 
FCart=.FALSE. 
574 
List=0 
575 
renum=.TRUE. 
576 
OptReac=.FALSE. 
577 
OptProd=.FALSE. 
578 
CalcEReac=.FALSE. 
579 
CalcEProd=.FALSE. 
580 
EReac=0.d0 
581 
EProd=0.d0 
582 
OptGeom=1 
583 
GeomFile="EMPTY" 
584 
AnaGeom=.FALSE. 
585 
AnaFile="EMPTY" 
586 
GplotFile="EMPTY" 
587 
Nb=0 
588 
NbCom=0 
589 
PathOnly=.FALSE. 
590 
Prog='EMPTY' 
591 
ProgExe='EMPTY' 
592 
Runmode='SERIAL' 
593 
Input='EMPTY' 
594 
Poscar="POSCAR" 
595 
AutoCart=.FALSE. 
596 
VMD=.TRUE. 
597 
WriteVASP=.FALSE. 
598 
DynMaxStep=.FALSE. 
599 
CalcName="EMPTY" 
600 
ISuffix="EMPTY" 
601 
OSuffix="EMPTY" 
602 
NGintMax=10 
603 
Align=.TRUE. 
604  
605 
! Number of point used to compute the length of the path, 
606 
! and to reparameterize the path 
607 
NMaxPtPath=1 
608 
FrozTol=1e4 
609  
610 
! Read the namelist 
611 
read (IOIN,path) 
612  
613 
debug=valid("Main").or.valid("Path") 
614  
615 
FCalcHess=valid("CalcHess") 
616 
PathName=AdjustL(PathName) 
617 
Coord=AdjustL(Coord) 
618 
CALL UpCase(coord) 
619  
620 
Step_method=AdjustL(Step_method) 
621 
CALL UpCase(Step_method) 
622  
623 
Prog=AdjustL(Prog) 
624 
CALL UpCase(prog) 
625  
626 
Input=AdjustL(Input) 
627 
CALL UpCase(input) 
628  
629 
Runmode=AdjustL(Runmode) 
630 
Call UpCase(Runmode) 
631 
IF (Runmode(1:4).EQ.'PARA') Runmode="PARA" 
632  
633 
if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN 
634 
WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'" 
635 
WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE" 
636 
CalcEReac=.FALSE. 
637 
CalcEProd=.FALSE. 
638 
END IF 
639  
640 
If (IReparamT<=0) IReparamT=IReparam 
641  
642 
! We take care of HesUpd flag 
643 
! this is needed because some users might still use it 
644 
if (HesUpd.NE."EMPTY") THEN 
645 
! We are pityless: 
646 
WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead" 
647 
WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead" 
648 
STOP 
649 
END IF 
650  
651 
IF (HUpdate=='EMPTY') THEN 
652 
IF (OptGeom>=1) THEN 
653 
HupDate="BFGS" 
654 
ELSE 
655 
HUpdate="MS" 
656 
END IF 
657 
END IF 
658  
659 
If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 
660 
COORD='CART' 
661 
END IF 
662  
663 
IF (COORD.EQ.'CART') THEN 
664 
Renum=.FALSE. 
665 
IF (OptGeom>0) THEN 
666 
IGeomRef=OptGeom 
667 
ELSE 
668 
IGeomRef=1 
669 
END IF 
670 
END IF 
671  
672 

673 
if (Prog.EQ.'TEST2') Prog="CHAMFRE" 
674 
if (Prog.EQ.'2D') Prog="TEST2D" 
675 
if (Prog.EQ.'TEST3') Prog="LEPS" 
676  
677  
678 
Select Case (Prog) 
679 
CASE ("EMPTY","G09","GAUSSIAN") 
680 
Prog='GAUSSIAN' 
681 
If (ProgExe=="EMPTY") ProgExe="g09" 
682 
if (CalcName=="EMPTY") CalcName="Path" 
683 
if (ISuffix=="EMPTY") ISuffix="com" 
684 
if (OSuffix=="EMPTY") OSuffix="log" 
685 
UnitProg="au" 
686 
ConvE=au2kcal 
687 
CASE ("G03") 
688 
Prog='GAUSSIAN' 
689 
If (ProgExe=="EMPTY") ProgExe="g03" 
690 
if (CalcName=="EMPTY") CalcName="Path" 
691 
if (ISuffix=="EMPTY") ISuffix="com" 
692 
if (OSuffix=="EMPTY") OSuffix="log" 
693 
UnitProg="au" 
694 
ConvE=au2kcal 
695 
CASE ("EXT") 
696 
If (ProgExe=="EMPTY") ProgExe="./Ext.exe" 
697 
if (CalcName=="EMPTY") CalcName="Ext" 
698 
if (ISuffix=="EMPTY") ISuffix="in" 
699 
if (OSuffix=="EMPTY") OSuffix="out" 
700 
UnitProg="au" 
701 
ConvE=au2kcal 
702 
CASE ('MOPAC') 
703 
If (ProgExe=="EMPTY") ProgExe="mopac" 
704 
if (CalcName=="EMPTY") CalcName="Path" 
705 
if (ISuffix=="EMPTY") ISuffix="mop" 
706 
if (OSuffix=="EMPTY") OSuffix="out" 
707 
UnitProg="au" 
708 
ConvE=au2kcal 
709 
CASE ("SIESTA") 
710 
If (ProgExe=="EMPTY") ProgExe="siesta" 
711 
if (CalcName=="EMPTY") CalcName="Path" 
712 
if (ISuffix=="EMPTY") ISuffix="fdf" 
713 
if (OSuffix=="EMPTY") OSuffix="out" 
714 
UnitProg="eV" 
715 
ConvE=eV2kcal 
716 
CASE ("VASP") 
717 
If (ProgExe=="EMPTY") ProgExe="VASP" 
718 
UnitProg="eV" 
719 
ConvE=eV2kcal 
720 
CASE ("TM","TURBOMOLE") 
721 
Prog="TURBOMOLE" 
722 
If (ProgExe=="EMPTY") ProgExe="jobex c 1 keep" 
723 
UnitProg="au" 
724 
ConvE=au2kcal 
725 
CASE ("TEST","CHAMFRE","TEST2D","LEPS") 
726 
ProgExe="" 
727 
UnitProg="au" 
728 
ConvE=au2kcal 
729 
CASE DEFAULT 
730 
WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP" 
731 
STOP 
732 
END SELECT 
733  
734 
if (Input.EQ.'EMPTY') THEN 
735 
Select Case (Prog) 
736 
CASE ("VASP") 
737 
Input=Prog 
738 
CASE ("CHAMFRE") 
739 
Input="VASP" 
740 
CASE DEFAULT 
741 
Input='XYZ' 
742 
END SELECT 
743 
WRITE(*,*) "Input has been set to the default: ",INPUT 
744 
END IF 
745  
746 
WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) 
747  
748 
IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN 
749 
WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode." 
750 
STOP 
751 
END IF 
752  
753 
if (OptGeom.GE.1) THEN 
754 
FPrintGeom=.FALSE. 
755 
If (GeomFile/='EMPTY') THEN 
756 
FPrintGeom=.TRUE. 
757 
END IF 
758 
IF (IGeomRef.LE.0) THEN 
759 
IGeomRef=OptGeom 
760 
ELSE 
761 
IF (IGeomRef.NE.OptGeom) THEN 
762 
WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
763 
WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
764  
765 
WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
766 
WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
767 
STOP 
768 
END IF 
769 
END IF 
770 
END IF 
771  
772 
ALLOCATE(Cart(Nat)) 
773 
Cart=0 
774  
775 
ALLOCATE(Frozen(Nat)) 
776 
Frozen=0 
777  
778 
IF (FCart) THEN 
779 
! We rewind the file to be sure to browse all of it to read the Cart namelist 
780 
REWIND(IOIN) 
781 
List=0 
782 
READ(IOIN,cartlist) 
783 
IF (COORD.NE.'CART') THEN 
784 
Cart=List(1:Nat) 
785 
if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart 
786 
ELSE 
787 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
788 
WRITE(*,*) "I will discard the cartlist" 
789 
Cart=0 
790 
END IF 
791 
END IF 
792  
793 
if (FFrozen) THEN 
794  
795 
REWIND(IOIN) 
796 
List=0 
797 
READ(IOIN,Frozenlist) 
798 
Frozen=List(1:Nat) 
799 
if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen 
800 
END IF 
801  
802 
If (FCart) Call Expand(Cart,NCart,Nat) 
803 
IF (FFrozen) Call Expand(Frozen,NFroz,Nat) 
804  
805 
if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart 
806 
if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen 
807  
808 
if (AnaGeom) THEN 
809 
REWIND(IOIN) 
810 
READ(IOIN,AnaList,IOSTAT=IoS) 
811 
IF (IOS==0) THEN ! we have a AnaList namelist 
812 
Call ReadAnaList 
813 
ELSE 
814 
! No AnaList namelist, we will print only the energy 
815 
FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)' 
816 
END IF 
817 
IF (AnaFile=="EMPTY") THEN 
818 
AnaFile=TRIM(PathName) // '.dat' 
819 
END IF 
820 
OPEN(IODat,File=AnaFile) 
821 
IF (GplotFile=="EMPTY") THEN 
822 
GplotFile=TRIM(PathName) // '.gplot' 
823 
END IF 
824 

825 
OPEN(IODat,File=AnaFile) 
826 
If (OptGeom<=0) THEN 
827 
OPEN(IOGplot,File=GPlotFile) 
828 
WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot persist" 
829 
WRITE(IoGplot,'(A)') " set pointsize 2" 
830 
WRITE(IoGplot,'(A)') " xcol=1" 
831 
WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3 
832 
END IF 
833 
END IF 
834  
835 

836  
837 
if (NMaxPtPath.EQ.1) then 
838 
NMaxPtPath=min(50*NGeomF,2000) 
839 
NMaxPtPath=max(20*NGeomF,NMaxPtPath) 
840 
end if 
841  
842 
! NMaxPtPath=10000 
843  
844 
! We ensure that FTan is in [0.:1.] 
845 
FTan=min(abs(FTan),1.d0) 
846  
847 
! PFL 2011 Mar 14 > 
848 
! Added some printing for analyses with Anapath 
849 
if (debug) THEN 
850 
WRITE(IOOUT,path) 
851 
ELSE 
852 
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC 
853 
! and PATHNAME 
854 
WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat 
855 
WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI 
856 
WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF 
857 
WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC 
858 
WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME) 
859 
END IF 
860  
861 
FTan=1.FTan 
862  
863 
!Thresholds... 
864 
if (SThresh.LE.0) SThresh=0.01d0 
865 
If (GThresh.LE.0) GThresh=SThresh/4. 
866 
SThreshRms=Sthresh*2./3. 
867 
GThreshRms=Gthresh*2./3. 
868  
869 
! First batch of allocations. Arrays that depend only on NGeomI, Nat 
870 
ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) 
871 
ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat)) 
872 
ALLOCATE(MassAt(NAt),OrderInv(Nat)) 
873  
874 
AtName="" 
875 
MassAt=1. 
876  
877 
! As we have used rewind many times, the position in the input file 
878 
! is blurry (at least !) we thus have to go back to the Geometries description 
879 
! (TO DO: use a better Input parser !) 
880  
881 
REWIND(IOIN) 
882 
! We search the '&path' namelist 
883 
READ(IOIN,'(A)',Iostat=Ios) Line 
884 
Fini=.FALSE. 
885 
DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0) 
886 
if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line) 
887 
READ(IOIN,'(A)',Iostat=Ios) Line 
888 
END DO 
889 
if (debug) WRITE(*,*) "Path found:",TRIM(LINE) 
890 
! We are on the &path line, we move to the end 
891 
Call NoComment(Line) 
892 
DO WHILE ((Ios==0).AND.InString(Line,"/")==0) 
893 
if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line) 
894 
READ(IOIN,'(A)',Iostat=Ios) Line 
895 
Call NoComment(Line) 
896 
Call noString(Line) 
897 
END DO 
898 
! We are at then end of &Path / namelist 
899 
! We scan for other Namelists 
900 
READ(IOIN,'(A)',Iostat=Ios) Line 
901 
if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line) 
902 
DO WHILE ((IoS==0).AND.(Index(Line,'&')/=0)) 
903 
Call NoComment(Line) 
904 
if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line) 
905 
Idx=InString(Line,'Analist') 
906 
DO WHILE ((Ios==0).AND.InString(Line,"/")==0) 
907 
if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line) 
908 
READ(IOIN,'(A)',Iostat=Ios) Line 
909 
Call NoComment(Line) 
910 
END DO 
911 
If (Idx/=0) THEN 
912 
! we have just read &Analist namelist, we have to skip the variables description 
913 
DO I=1,Nb 
914 
READ(IOIN,'(A)',Iostat=Ios) Line 
915 
if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE) 
916 
END DO 
917 
END IF 
918 
READ(IOIN,'(A)',Iostat=Ios) Line 
919 
if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE) 
920 
END DO 
921 
BACKSPACE(IOIN) 
922  
923 
! We read the initial geometries 
924 
Call Read_Geom(input) 
925  
926 
! We convert atome names into atom mass number 
927 
DO I=1,NAt 
928 
! WRITE(*,*) 'I,AtName',I,AtName(I) 
929 
Atome(I)=ConvertNumAt(AtName(I)) 
930 
END DO 
931  
932 
If (AnaGeom) THEN 
933 
If (NbVar>0) THEN 
934 
! We print the description of the varialbes in AnaFile 
935 
Call PrintAnaList(IoDat) 
936 
If (OptGeom<=0) THEN 
937 
WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")" 
938 
ELSE 
939 
WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")" 
940 
END IF 
941 
ELSE 
942 
If (OptGeom<=0) THEN 
943 
WRITE(IoDat,'(A)') "# s Erel(kcal/mol) E (" // TRIM(UnitProg) //")" 
944 
ELSE 
945 
WRITE(IoDat,'(A)') "# Iter Erel(kcal/mol) E (" // TRIM(UnitProg) //")" 
946 
END IF 
947 
END IF 
948 
END IF 
949  
950 
! PFL 23 Sep 2008 > 
951 
! We check that frozen atoms are really frozen, ie that their coordinates do not change 
952 
! between the first geometry and the others. 
953 
IF (NFroz.GT.0) THEN 
954 
ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat)) 
955 
ALLOCATE(x2(nat),y2(nat),z2(nat)) 
956 
ListAtFroz=.FALSE. 
957  
958 
x1(1:Nat)=XyzgeomI(1,1,1:Nat) 
959 
y1(1:Nat)=XyzgeomI(1,2,1:Nat) 
960 
z1(1:Nat)=XyzgeomI(1,3,1:Nat) 
961  
962 
SELECT CASE (NFroz) 
963 
! We have Frozen atoms 
964 
! PFL 28 Oct 2008 > 
965 
! It might happen that the geometries are translated/rotated. 
966 
! So if we have 3 (or more) frozen atoms, we reorient the geometries, based on the first 
967 
! three frozen atoms. 
968  
969 
CASE (3:) 
970 
DO I=1,NFroz 
971 
ListAtFroz(Frozen(I))=.TRUE. 
972 
END DO 
973 
CASE (2) 
974 
IAt1=Frozen(1) 
975 
IAt2=Frozen(2) 
976 
ListAtFroz(Iat1)=.TRUE. 
977 
ListAtFroz(Iat2)=.TRUE. 
978 
x2(1)=x1(Iat2)x1(Iat1) 
979 
y2(1)=y1(Iat2)y1(Iat1) 
980 
z2(1)=z1(Iat2)z1(Iat1) 
981 
Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) 
982 
Dist=1. 
983 
! We scan all atoms to find one that is not aligned with Iat1Iat2 
984 
DO I=1,Nat 
985 
if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle 
986 
x2(2)=x1(Iat2)x1(I) 
987 
y2(2)=y1(Iat2)y1(I) 
988 
z2(2)=z1(Iat2)z1(I) 
989 
Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2) 
990 
ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2) 
991 
if (abs(ca)<0.9) THEN 
992 
IF (Norm2>Dist) THEN 
993 
Iat3=I 
994 
Dist=Norm2 
995 
END IF 
996 
END IF 
997 
END DO 
998 
ListAtFroz(Iat3)=.TRUE. 
999 
if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3 
1000 
CASE (1) 
1001 
IAt1=Frozen(1) 
1002 
ListAtFroz(Iat1)=.TRUE. 
1003 
Dist=1. 
1004 
! We scan all atoms to find one that is further from At1 
1005 
DO I=1,Nat 
1006 
if (I.EQ.Iat1) Cycle 
1007 
x2(1)=x1(Iat1)x1(I) 
1008 
y2(1)=y1(Iat1)y1(I) 
1009 
z2(1)=z1(Iat1)z1(I) 
1010 
Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) 
1011 
IF (Norm1>Dist) THEN 
1012 
Iat2=I 
1013 
Dist=Norm1 
1014 
END IF 
1015 
END DO 
1016 
ListAtFroz(Iat2)=.TRUE. 
1017 
if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2 
1018 
x2(1)=x1(Iat2)x1(Iat1) 
1019 
y2(1)=y1(Iat2)y1(Iat1) 
1020 
z2(1)=z1(Iat2)z1(Iat1) 
1021 
Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) 
1022 
Dist=1. 
1023 
! We scan all atoms to find one that is not aligned with Iat1Iat2 
1024 
Iat3=1 
1025 
DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat)) 
1026 
Iat3=Iat3+1 
1027 
END DO 
1028 
DO I=1,Nat 
1029 
if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle 
1030 
x2(2)=x1(Iat2)x1(I) 
1031 
y2(2)=y1(Iat2)y1(I) 
1032 
z2(2)=z1(Iat2)z1(I) 
1033 
Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2) 
1034 
ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2) 
1035 
if (abs(ca)<0.9) THEN 
1036 
IF (Norm2>Dist) THEN 
1037 
Iat3=I 
1038 
Dist=Norm2 
1039 
END IF 
1040 
END IF 
1041 
END DO 
1042 
ListAtFroz(Iat3)=.TRUE. 
1043 
if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3 
1044 
END SELECT 
1045  
1046 
DO I=2,NGeomI 
1047 
! First, we align the Ith geometry on the first one, using only 
1048 
! the positions of the "frozen" atoms. (see above: it might be that 
1049 
! ListAtFroz contains non frozen atoms used to align the geometries 
1050 
x2(1:Nat)=XyzgeomI(I,1,1:Nat) 
1051 
y2(1:Nat)=XyzgeomI(I,2,1:Nat) 
1052 
z2(1:Nat)=XyzgeomI(I,3,1:Nat) 
1053 
Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot) 
1054  
1055 
FChkFrozen=.FALSE. 
1056 
DO J=1,NFroz 
1057 
IAt=Frozen(J) 
1058 
IF ((abs(X1(Iat)X2(Iat)).GT.FrozTol).OR. & 
1059 
(abs(y1(Iat)y2(Iat)).GT.FrozTol).OR. & 
1060 
(abs(z1(Iat)z2(Iat)).GT.FrozTol)) THEN 
1061 
IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****" 
1062 
FChkFrozen=.TRUE. 
1063 
WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change" 
1064 
WRITE(*,*) "from geometry 1 to geometry ",I 
1065 
END IF 
1066 
END DO 
1067 
END DO 
1068 
IF (FChkFrozen) THEN 
1069 
WRITE(*,*) "Some frozen atoms are not frozen... check this and play again" 
1070 
STOP 
1071 
END IF 
1072 
END IF 
1073  
1074  
1075 
N3at=3*Nat 
1076 
NFree=N3at6 ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES. 
1077  
1078  
1079 
Call ReadInput 
1080  
1081 
IF (COORD=="MIXED") Call TestCart(AutoCart) 
1082  
1083 
SELECT CASE(NFroz) 
1084 
CASE (2) 
1085 
IntFroz=1 
1086 
CASE (3) 
1087 
IntFroz=3 
1088 
CASE (4:) 
1089 
IntFroz=(NFroz2)*3 !(NFroz3)*3+3 
1090 
END SELECT 
1091  
1092 
if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "& 
1093 
,TRIM(COORD),Ncart,NCoord,N3at,NFroz 
1094 
if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart 
1095  
1096 
if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN 
1097 
WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT." 
1098 
COORD="ZMAT" 
1099 
END IF 
1100  
1101 
IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz 
1102  
1103 
SELECT CASE(COORD) 
1104 
CASE ('ZMAT') 
1105 
NCoord=Nfree 
1106 
ALLOCATE(dzdc(3,nat,3,nat)) 
1107 
CASE ('MIXED') 
1108 
SELECT CASE (NCart+NFroz) 
1109 
CASE (1) 
1110 
NCoord=N3At3 
1111 
CASE (2) 
1112 
NCoord=N3At1 
1113 
CASE (3:) 
1114 
NCoord=N3At 
1115 
END SELECT 
1116 
ALLOCATE(dzdc(3,nat,3,nat)) 
1117 
CASE ('BAKER') 
1118 
Symmetry_elimination=0 
1119 
NCoord=3*Nat6Symmetry_elimination !NFree = 3*Nat6 
1120 
ALLOCATE(BMat_BakerT(3*Nat,NCoord)) 
1121 
ALLOCATE(BTransInv(NCoord,3*Nat)) 
1122 
ALLOCATE(BTransInv_local(NCoord,3*Nat)) 
1123 
CASE ('HYBRID') 
1124 
NCoord=N3at 
1125 
CASE ('CART') 
1126 
NCoord=N3at 
1127 
END SELECT 
1128  
1129 
if (debug) WRITE(*,*) "DBG Path, L826, Coord, NCart,NCoord, N3At",TRIM(COORD),Ncart, NCoord,N3at 
1130  
1131 
ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord)) 
1132 
ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat)) 
1133 
ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord)) 
1134 
ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in Path_module.f90 
1135 
ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord)) 
1136 
ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord)) 
1137 
ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in Path_module.f90 
1138 
ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF)) 
1139  
1140 
ZeroVec=0.d0 
1141 
DO I =1, NGeomF 
1142 
IntTangent(I,:)=0.d0 
1143 
END DO 
1144  
1145 
if (debug) THEN 
1146 
Print *, 'L886, IntTangent(1,:)=' 
1147 
Print *, IntTangent(1,:) 
1148 
END IF 
1149  
1150 
if (.NOT.OptReac) Energies(1)=Ereac 
1151 
if (.NOT.OptProd) Energies(NGeomF)=EProd 
1152 
MaxStep=SMax 
1153  
1154 
! We initialize the mass array 
1155 
if (MW) THEN 
1156 
WRITE(*,*) 'Working in MW coordinates' 
1157 
DO I=1,Nat 
1158 
massAt(I)=Mass(Atome(I)) 
1159 
END DO 
1160 
ELSE 
1161 
DO I=1,Nat 
1162 
MassAt(I)=1.d0 
1163 
END DO 
1164 
END IF 
1165  
1166 
WRITE(*,*) "Prog=",TRIM(PROG) 
1167  
1168 
ALLOCATE(FrozAtoms(Nat)) 
1169  
1170 
! IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz 
1171 
! IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen 
1172 
! IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms 
1173 
IF (NFroz.EQ.0) THEN 
1174 
FrozAtoms=.TRUE. 
1175 
ELSE 
1176 
I=1 
1177 
NFr=0 
1178 
FrozAtoms=.False. 
1179 
DO WHILE (Frozen(I).NE.0) 
1180 
IF (Frozen(I).LT.0) THEN 
1181 
DO J=Frozen(I1),abs(Frozen(I)) 
1182 
IF (.NOT.FrozAtoms(J)) THEN 
1183 
NFr=NFr+1 
1184 
FrozAtoms(J)=.TRUE. 
1185 
END IF 
1186 
END DO 
1187 
ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN 
1188 
FrozAtoms(Frozen(I))=.TRUE. 
1189 
NFr=NFr+1 
1190 
END IF 
1191 
I=I+1 
1192 
END DO 
1193 
IF (NFr.NE.NFroz) THEN 
1194 
WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :( ",NFroz,NFr 
1195 
STOP 
1196 
END IF 
1197 
END IF 
1198  
1199 
! if (debug) THEN 
1200 
!Print *, 'L968, IntTangent(1,:)=' 
1201 
!Print *, IntTangent(1,:) 
1202 
! END IF 
1203  
1204 
! We have read everything we needed... time to work :) 
1205 
IOpt=0 
1206 
FirstTimePathCreate = .TRUE. ! Initialized. 
1207 
Call PathCreate(IOpt) ! IntTangent is also computed in PathCreate. 
1208 
FirstTimePathCreate = .FALSE. 
1209  
1210 
if (debug) THEN 
1211 
Print *, 'L980, IntTangent(1,:)=' 
1212 
Print *, IntTangent(1,:) 
1213 
END IF 
1214  
1215 
! PFL 30 Mar 2008 > 
1216 
! The following is now allocated in Calc_Baker. This is more logical to me 
1217 
! IF (COORD .EQ. "BAKER") Then 
1218 
! ALLOCATE(XprimitiveF(NgeomF,NbCoords)) 
1219 
! ALLOCATE(UMatF(NGeomF,NbCoords,NCoord)) 
1220 
! ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat)) 
1221 
! END IF 
1222 
! < PFL 30 mar 2008 
1223  
1224 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1225 
! 
1226 
! For Debugging purposes 
1227 
! 
1228 
! OptGeom is the index of the geometry which is to be optimized. 
1229 
IF (Optgeom.GT.0) THEN 
1230 
Flag_Opt_Geom=.TRUE. 
1231 
SELECT CASE(Coord) 
1232 
CASE ('CART','HYBRID') 
1233 
!!! CAUTION : PFL 29.AUG.2008 > 
1234 
! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3) 
1235 
! This should be made consistent !!!!!!!!!!!!!!! 
1236 
GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
1237 
! GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1238 
!!! < CAUTION : PFL 29.AUG.2008 
1239 
CASE ('ZMAT','MIXED') 
1240 
!GeomTmp=IntCoordF(OptGeom,:) 
1241 
GeomTmp=IntCoordI(OptGeom,:) 
1242 
CASE ('BAKER') 
1243 
!GeomTmp=IntCoordF(OptGeom,:) 
1244 
GeomTmp=IntCoordI(OptGeom,:) 
1245 
CASE DEFAULT 
1246 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP" 
1247 
STOP 
1248 
END SELECT 
1249  
1250 
! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker 
1251 
Flag_Opt_Geom = .TRUE. 
1252 
Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method) 
1253  
1254 
WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom 
1255 
WRITE(*,*) "Before PrintGeom, L1059, Path.f90" 
1256 
Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat) 
1257 
STOP 
1258 
END IF ! IF (Optgeom.GT.0) THEN 
1259  
1260 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1261 
! End of GeomOpt 
1262 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1263  
1264 
IF (PathOnly) THEN 
1265 
WRITE(*,*) "PathOnly=.T. , Stop here" 
1266 
Call Write_path(1) 
1267 
if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar) 
1268 
STOP 
1269 
END IF 
1270  
1271 
if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) 
1272  
1273 
DEALLOCATE(XyzGeomI,IntCoordI) 
1274 
NGeomI=NGeomF 
1275 
ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord)) 
1276  
1277 
IF (DEBUG) WRITE(*,*) ' Back from PathCreate, L965, Path.f90' 
1278 
! we print this path, which is the first one :) called Iteration0 
1279 
! we now have to calculate the gradient for each point (and the energy 
1280 
! if we work at 0K) 
1281  
1282 
if (DEBUG.AND.(COORD.EQ."BAKER")) THEN 
1283 
DO IGeom=1,NGeomF 
1284 
WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom 
1285 
DO I=1,3*Nat6 
1286 
WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) 
1287 
END DO 
1288 
END DO 
1289 
END IF 
1290  
1291 
! To calculate B^T^1 for all extrapolated final geometries: 
1292 
! PFL 18 Jan 2008 > 
1293 
! Added a test for COORD=BAKER before the call to calc_baker_allgeomF 
1294 
IF (COORD.EQ."BAKER") Call Calc_baker_allGeomF() 
1295 
! < PFL 18 Jan 2008 
1296  
1297 
if (DEBUG.AND.(COORD.EQ."BAKER")) THEN 
1298 
DO IGeom=1,NGeomF 
1299 
WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom 
1300 
DO I=1,3*Nat6 
1301 
WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) 
1302 
END DO 
1303 
END DO 
1304 
END IF 
1305  
1306 
IOpt=0 
1307 
Call EgradPath(IOpt) 
1308  
1309 
if (debug) WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt 
1310 
Call Write_path(IOpt) 
1311 
! Shall we analyze the geometries ? 
1312 
IF (AnaGeom) Call AnaPath(Iopt,IoDat) 
1313  
1314 
if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) 
1315  
1316 
! We have the geometries, the first gradients... we will generate the first hessian matrices :) 
1317 
! ... v3.71 Only if IniHUp=.TRUE. which is not the default. In this version, we take the hessian 
1318 
! of end points as Unity matrices, and we update them along the path. 
1319  
1320 
! By default, Hess=Id 
1321 
Hess=0. 
1322 
IF(Coord .EQ. "ZMAT") Then 
1323 
! We use the fact that we know that approximate good hessian values 
1324 
! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles 
1325 
DO IGeom=1,NGeomF 
1326 
Hess(IGeom,1,1)=0.5d0 
1327 
Hess(IGeom,2,2)=0.5d0 
1328 
Hess(IGeom,3,3)=0.1d0 
1329 
DO J=1,Nat3 
1330 
Hess(IGeom,3*J+1,3*J+1)=0.5d0 
1331 
Hess(IGeom,3*J+2,3*J+2)=0.1d0 
1332 
Hess(IGeom,3*J+3,3*J+3)=0.02d0 
1333 
END DO 
1334 
IF (HInv) THEN 
1335 
DO I=1,NCoord 
1336 
Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I) 
1337 
END DO 
1338 
END IF 
1339 
END DO 
1340 
ELSE 
1341 
DO I=1,NCoord 
1342 
DO J=1,NGeomF 
1343 
Hess(J,I,I)=1.d0 
1344 
END DO 
1345 
END DO 
1346 
END IF 
1347  
1348 
IF (COORD .EQ. "BAKER") THEN 
1349 
! UMat(NPrim,NCoord) 
1350 
ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat6)) 
1351 
Hprim=0.d0 
1352 
ScanCoord=>Coordinate 
1353 
I=0 
1354 
DO WHILE (Associated(ScanCoord%next)) 
1355 
I=I+1 
1356 
SELECT CASE (ScanCoord%Type) 
1357 
CASE ('BOND') 
1358 
!DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF?? 
1359 
Hprim(I,I) = 0.5d0 
1360 
!END DO 
1361 
CASE ('ANGLE') 
1362 
Hprim(I,I) = 0.2d0 
1363 
CASE ('DIHEDRAL') 
1364 
Hprim(I,I) = 0.1d0 
1365 
END SELECT 
1366 
ScanCoord => ScanCoord%next 
1367 
END DO 
1368  
1369 
! HprimU = H^prim U: 
1370 
HprimU=0.d0 
1371 
DO I=1, 3*Nat6 
1372 
DO J=1,NPrim 
1373 
HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I) 
1374 
END DO 
1375 
END DO 
1376  
1377 
! Hess = U^T HprimU = U^T H^prim U: 
1378 
Hess=0.d0 
1379 
DO I=1, 3*Nat6 
1380 
DO J=1,NPrim 
1381 
! UMat^T is needed below. 
1382 
Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I) 
1383 
END DO 
1384 
END DO 
1385 
DO K=2,NGeomF 
1386 
Hess(K,:,:)=Hess(1,:,:) 
1387 
END DO 
1388 
DEALLOCATE(Hprim,HprimU) 
1389 
end if ! ends IF COORD=BAKER 
1390 
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord)) 
1391 
ALLOCATE(HessTmp(NCoord,NCoord)) 
1392  
1393 
IF (IniHUp) THEN 
1394 
IF (FCalcHess) THEN 
1395 
! The numerical calculation of the Hessian has been switched on 
1396 
DO IGeom=1,NGeomF 
1397 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1398 
GeomTmp=IntCoordF(IGeom,:) 
1399 
ELSE 
1400 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1401 
END IF 
1402 
Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt) 
1403 
END DO 
1404 
ELSE 
1405 
! We generate a forward hessian and a backward one and we mix them. 
1406 
! First, the forward way: 
1407 
DO IGeom=2,NGeomF1 
1408 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1409 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom1,:) 
1410 
ELSE 
1411 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom1,:,:),(/3*Nat/)) 
1412 
END IF 
1413 
GradTmp=Grad(IGeom1,:)Grad(IGeom,:) 
1414 
HessTmp=Hess(IGeom1,:,:) 
1415 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
1416 
Hess(IGeom,:,:)=HessTmp 
1417 
END DO 
1418  
1419 
! Now, the backward way: 
1420 
HessTmp=Hess(NGeomF,:,:) 
1421 
DO IGeom=NGeomF1,2,1 
1422 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1423 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom+1,:) 
1424 
ELSE 
1425 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) 
1426 
END IF 
1427 
GradTmp=Grad(IGeom,:)Grad(IGeom+1,:) 
1428 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
1429 
alpha=0.5d0*Pi*(IGeom1.)/(NGeomF1.) 
1430 
ca=cos(alpha) 
1431 
sa=sin(alpha) 
1432 
Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:) 
1433 
END DO 
1434 
END IF ! matches IF (FCalcHess) 
1435 
END IF ! matches IF (IniHUp) THEN 
1436  
1437 
! For debug purposes, we diagonalize the Hessian matrices 
1438 
IF (debug) THEN 
1439 
!WRITE(*,*) "DBG Main, L1104,  EigenSystem for Hessian matrices" 
1440 
DO I=1,NGeomF 
1441 
WRITE(*,*) "DBG Main  Point ",I 
1442 
Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord) 
1443 
DO J=1,NCoord 
1444 
WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord)) 
1445 
END DO 
1446 
END DO 
1447 
END IF ! matches IF (debug) THEN 
1448  
1449 
! we have the hessian matrices and the gradients, we can start the optimizations ! 
1450 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1451 
! 
1452 
! Main LOOP : Optimization of the Path 
1453 
! 
1454 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1455 
if (debug) Print *, 'NGeomF=', NGeomF 
1456 
allocation_flag = .TRUE. 
1457  
1458 
Fini=.FALSE. 
1459  
1460 
DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini)) 
1461 
if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt 
1462 
IOpt=IOpt+1 
1463  
1464 
SELECT CASE (COORD) 
1465 
CASE ('ZMAT','MIXED','BAKER') 
1466 
GeomOld_all=IntCoordF 
1467 
CASE ('CART','HYBRID') 
1468 
GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/)) 
1469 
CASE DEFAULT 
1470 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP' 
1471 
STOP 
1472 
END SELECT 
1473  
1474 
IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. & 
1475 
valid("printtangent")) THEN 
1476 
WRITE(*,*) "Visualization of tangents" 
1477 
Idx=min(12,NCoord) 
1478 
DO I=1,NGeomF 
1479 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)0.5d0*IntTangent(I,1:Idx) 
1480 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx) 
1481 
WRITE(*,*) 
1482 
END DO 
1483 
WRITE(*,*) "END of tangents" 
1484 
END IF 
1485  
1486 
IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. & 
1487 
valid("printgrad")) THEN 
1488 
WRITE(*,*) "Visualization of gradients" 
1489 
Idx=min(12,NCoord) 
1490 
DO I=1,NGeomF 
1491 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)1.d0*Grad(I,1:Idx) 
1492 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx) 
1493 
WRITe(*,*) 
1494 
END DO 
1495 
WRITE(*,*) "END of gradients" 
1496 
END IF 
1497  
1498 
Fini=.TRUE. 
1499 
IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE. 
1500 
IGeom=1 
1501 
if (debug) WRITE(*,*) '**************************************' 
1502 
if (debug) WRITE(*,*) 'Optimizing reactant' 
1503 
if (debug) WRITE(*,*) '**************************************' 
1504 
SELECT CASE (COORD) 
1505 
CASE ('ZMAT','MIXED','BAKER') 
1506 
SELECT CASE (Step_method) 
1507 
CASE ('RFO') 
1508 
!GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT). 
1509 
Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec) 
1510 
Print *, GradTmp 
1511 
CASE ('GDIIS') 
1512 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).: 
1513 
Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,& 
1514 
Hess(1,:,:),NCoord,allocation_flag,ZeroVec) 
1515 
Print *, 'OptReac, Steps from GDIIS, GeomNew  IntCoordF(1,:)=' 
1516 
Print *, GradTmp 
1517 
CASE ('GEDIIS') 
1518 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).: 
1519 
! Energies are updated in EgradPath.f90 
1520 
Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),& 
1521 
NCoord,allocation_flag,ZeroVec) 
1522 
!Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag) 
1523 
!GradTmp = GeomTmp  IntCoordF(1,:) 
1524 
!Print *, 'Old Geometry:' 
1525 
!Print *, IntCoordF(1,:) 
1526 
Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp  IntCoordF(1,:)=' 
1527 
Print *, GradTmp 
1528 
!Print *, 'GeomTmp:' 
1529 
!Print *, GeomTmp 
1530 
GeomTmp(:) = GradTmp(:) + IntCoordF(1,:) 
1531 
CASE DEFAULT 
1532 
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop" 
1533 
STOP 
1534 
END SELECT 
1535 
CASE ('HYBRID','CART') 
1536 
Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec) 
1537 
CASE DEFAULT 
1538 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP' 
1539 
STOP 
1540 
END SELECT 
1541  
1542 
IF (debug) THEN 
1543 
WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1' 
1544 
Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
1545 
END IF 
1546  
1547 
! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax 
1548 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 
1549 
FactStep=min(1.d0,MaxStep(Igeom)/NormStep) 
1550 
Fini=(Fini.AND.(NormStep.LE.SThresh)) 
1551 
OptReac=(NormStep.GT.SThresh) 
1552 
IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac 
1553 
NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/)))) 
1554 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
1555 
OptReac=(OptReac.OR.(NormGrad.GT.GThresh)) 
1556 
!Print *, 'Grad(1,:):' 
1557 
!Print *, Grad(1,:) 
1558 
IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac 
1559  
1560 
GradTmp=GradTmp*FactStep 
1561  
1562 
! we take care that frozen atoms do not move. 
1563 
IF (NFroz.NE.0) THEN 
1564 
SELECT CASE (COORD) 
1565 
CASE ('ZMAT','MIXED') 
1566 
GradTmp(1:IntFroz)=0.d0 
1567 
CASE ('CART','HYBRID') 
1568 
DO I=1,Nat 
1569 
IF (FrozAtoms(I)) THEN 
1570 
Iat=Order(I) 
1571 
GradTmp(3*Iat2:3*Iat)=0.d0 
1572 
END IF 
1573 
END DO 
1574 
CASE('BAKER') 
1575 
GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..) 
1576 
CASE DEFAULT 
1577 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP' 
1578 
STOP 
1579 
END SELECT 
1580 
END IF ! matches IF (NFroz.NE.0) THEN 
1581  
1582 
IF (debug) THEN 
1583 
!WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1' 
1584 
!Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
1585 
END IF 
1586  
1587 
IF (debug) THEN 
1588 
!WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1' 
1589 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1590 
WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree)) 
1591 
ELSE 
1592 
DO Iat=1,Nat 
1593 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1594 
XyzGeomF(IGeom,1:3,Iat) 
1595 
END DO 
1596 
END IF 
1597 
END IF ! matches IF (debug) THEN 
1598  
1599 
SELECT CASE (COORD) 
1600 
CASE ('ZMAT','MIXED','BAKER') 
1601 
IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1 
1602 
CASE ('HYBRID','CART') 
1603 
XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
1604 
CASE DEFAULT 
1605 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP' 
1606 
STOP 
1607 
END SELECT 
1608  
1609 
IF (debug) THEN 
1610 
WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1' 
1611 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1612 
WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree)) 
1613 
ELSE 
1614 
DO Iat=1,Nat 
1615  
1616 
!!!!!!!!!! There was an OrderInv here... should I put it back ? 
1617 
! It was with indzmat. IndZmat cannot appear here... 
1618 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1619 
XyzGeomF(IGeom,1:3,Iat) 
1620 
END DO 
1621 
END IF 
1622 
END IF ! matches IF (debug) THEN 
1623  
1624 
IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized" 
1625 
ELSE ! Optreac 
1626 
IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini 
1627 
END IF ! matches IF (OptReac) THEN 
1628  
1629 
IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE. 
1630 
IGeom=NGeomF 
1631 
if (debug) WRITE(*,*) '******************' 
1632 
if (debug) WRITE(*,*) 'Optimizing product' 
1633 
if (debug) WRITE(*,*) '******************' 
1634 
SELECT CASE (COORD) 
1635 
CASE ('ZMAT','MIXED','BAKER') 
1636 
Print *, 'Optimizing product' 
1637 
SELECT CASE (Step_method) 
1638 
CASE ('RFO') 
1639 
!GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT). 
1640 
Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec) 
1641 
Print *, GradTmp 
1642 
CASE ('GDIIS') 
1643 
! GradTmp becomes "step" below and has INTENT(OUT): 
1644 
Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),& 
1645 
HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec) 
1646 
Print *, 'OptProd, Steps from GDIIS, GeomNew  IntCoordF(NGeomF,:)=' 
1647 
Print *, GradTmp 
1648 
CASE ('GEDIIS') 
1649 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT): 
1650 
Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),& 
1651 
Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec) 
1652 
Print *, 'OptProd, Steps from GEDIIS, GeomNew  IntCoordF(NGeomF,:)=' 
1653 
Print *, GradTmp 
1654 
CASE DEFAULT 
1655 
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop" 
1656 
STOP 
1657 
END SELECT 
1658 
CASE ('HYBRID','CART') 
1659 
Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec) 
1660 
CASE DEFAULT 
1661 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP' 
1662 
STOP 
1663 
END SELECT 
1664  
1665 
! we have the Newton+RFO step, we will check that the displacement is not greater than Smax 
1666 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 
1667 
FactStep=min(1.d0,MaxStep(IGeom)/NormStep) 
1668 
Fini=(Fini.AND.(NormStep.LE.SThresh)) 
1669 
OptProd=.NOT.(NormStep.LE.SThresh) 
1670 
NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/)))) 
1671 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
1672 
OptProd=(OptProd.OR.(NormGrad.GT.GThresh)) 
1673 
IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd 
1674  
1675 
GradTmp=GradTmp*FactStep 
1676  
1677 
! we take care that frozen atoms do not move 
1678 
IF (NFroz.NE.0) THEN 
1679 
SELECT CASE (COORD) 
1680 
CASE ('ZMAT','MIXED','BAKER') 
1681 
GradTmp(1:IntFroz)=0.d0 
1682 
CASE ('CART','HYBRID') 
1683 
DO I=1,Nat 
1684 
IF (FrozAtoms(I)) THEN 
1685 
Iat=Order(I) 
1686 
GradTmp(3*Iat2:3*Iat)=0.d0 
1687 
END IF 
1688 
END DO 
1689 
CASE DEFAULT 
1690 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP' 
1691 
STOP 
1692 
END SELECT 
1693 
END IF ! matches IF (NFroz.NE.0) THEN 
1694  
1695 
IF (debug) THEN 
1696 
WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt 
1697 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1698 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1699 
ELSE 
1700 
DO Iat=1,Nat 
1701 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1702 
XyzGeomF(IGeom,1:3,Iat) 
1703 
END DO 
1704 
END IF 
1705 
END IF 
1706  
1707 
SELECT CASE (COORD) 
1708 
CASE ('ZMAT','MIXED','BAKER') 
1709 
IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF. 
1710 
CASE ('HYBRID','CART') 
1711 
XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
1712 
CASE DEFAULT 
1713 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP' 
1714 
STOP 
1715 
END SELECT 
1716  
1717 
IF (debug) THEN 
1718 
WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt 
1719 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1720 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1721 
ELSE 
1722 
DO Iat=1,Nat 
1723 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1724 
XyzGeomF(IGeom,1:3,Iat) 
1725 
END DO 
1726 
END IF 
1727 
END IF 
1728 
ELSE ! Optprod 
1729 
IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini 
1730 
END IF !matches IF (OptProd) THEN 
1731  
1732 
IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini 
1733  
1734 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1735 
! 
1736 
! Optimizing other geometries 
1737 
! 
1738 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1739  
1740  
1741  
1742 
DO IGeom=2,NGeomF1 ! matches at L1556 
1743 
if (debug) WRITE(*,*) '****************************' 
1744 
if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom 
1745 
if (debug) WRITE(*,*) '****************************' 
1746  
1747 
SELECT CASE (COORD) 
1748 
CASE ('ZMAT','BAKER','MIXED') 
1749 
GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/)) 
1750 
! PFL 6 Apr 2008 > 
1751 
! Special case, if FTan=0. we do not consider Tangent at all. 
1752 
! To DO: add the full treatment of FTan 
1753 
if (FTan==0.) GradTmp2=ZeroVec 
1754 
! < PFL 6 Apr 2008 
1755 
if (debug) THEN 
1756 
Print *, 'L1500, IntTangent(',IGeom,',:)=' 
1757 
Print *, IntTangent(IGeom,:) 
1758 
END IF 
1759 
!GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT). 
1760 
SELECT CASE (Step_method) 
1761 
CASE ('RFO') 
1762 
Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2) 
1763 
CASE ('GDIIS') 
1764 
!The following has no effect at IOpt==1 
1765 
!Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom 
1766 
Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,& 
1767 
Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2) 
1768 
Print *, 'All others, Steps from GDIIS, GeomNew  IntCoordF(IGeom,:)=' 
1769 
Print *, GradTmp 
1770 
CASE ('GEDIIS') 
1771 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT): 
1772 
Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),& 
1773 
Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2) 
1774 
Print *, 'All others, Steps from GEDIIS, GeomNew  IntCoordF(IGeom,:)=' 
1775 
Print *, GradTmp 
1776 
CASE DEFAULT 
1777 
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop" 
1778 
STOP 
1779 
END SELECT 
1780 
CASE ('HYBRID','CART') 
1781 
! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat 
1782 
! so we change it 
1783 
DO I=1,Nat 
1784 
DO J=1,3 
1785 
GradTmp2(J+(I1)*3)=XyzTangent(IGeom,(J1)*Nat+I) 
1786 
END DO 
1787 
END DO 
1788 
! GradTmp2=XyzTangent(IGeom,:) 
1789 
! PFL 6 Apr 2008 > 
1790 
! Special case, if FTan=0. we do not consider Tangent at all. 
1791 
! To DO: add the full treatment of FTan 
1792 
if (FTan==0.) GradTmp2=ZeroVec 
1793 
! < PFL 6 Apr 2008 
1794  
1795 
Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2) 
1796 
CASE DEFAULT 
1797 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP' 
1798 
STOP 
1799 
END SELECT 
1800  
1801 
! we take care that frozen atoms do not move 
1802 
IF (NFroz.NE.0) THEN 
1803 
SELECT CASE (COORD) 
1804 
CASE ('ZMAT','MIXED','BAKER') 
1805 
IF (debug) THEN 
1806 
WRITE(*,*) 'Step computed. Coord=',Coord 
1807 
WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz 
1808 
END IF 
1809 
GradTmp(1:IntFroz)=0.d0 
1810 
CASE ('CART','HYBRID') 
1811 
if (debug) THEN 
1812 
WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing" 
1813 
DO I=1,Nat 
1814 
WRITe(*,*) I,GradTmp(3*I2:3*I) 
1815 
END DO 
1816 
END IF 
1817 
DO I=1,Nat 
1818 
IF (FrozAtoms(I)) THEN 
1819 
if (debug) THEN 
1820 
write(*,*) 'Step Computed. Coord=',Coord 
1821 
WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I) 
1822 
END IF 
1823 
Iat=Order(I) 
1824 
GradTmp(3*Iat2:3*Iat)=0.d0 
1825 
END IF 
1826 
END DO 
1827  
1828 
if (debug) THEN 
1829 
WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing" 
1830 
DO I=1,Nat 
1831 
WRITe(*,*) I,GradTmp(3*I2:3*I) 
1832 
END DO 
1833 
END IF 
1834  
1835 
CASE DEFAULT 
1836 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP' 
1837 
STOP 
1838 
END SELECT 
1839 
END IF !matches IF (NFroz.NE.0) THEN 
1840  
1841 
IF (debug) THEN 
1842 
SELECT CASE (COORD) 
1843 
CASE ('ZMAT','MIXED','BAKER') 
1844 
WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree)) 
1845 
WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree)) 
1846 
CASE ('CART','HYBRID') 
1847 
WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree)) 
1848 
WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree)) 
1849 
CASE DEFAULT 
1850 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP' 
1851 
STOP 
1852 
END SELECT 
1853 
END IF ! matches if (debug) THEN 
1854  
1855 
! we check that the maximum displacement is not greater than Smax 
1856 
If (debug) WRITE(*,*) "Fini before test:",Fini 
1857 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 
1858 
FactStep=min(1.d0,maxStep(Igeom)/NormStep) 
1859 
Fini=(Fini.AND.(NormStep.LE.SThresh)) 
1860 
IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini 
1861  
1862 
GradTmp=GraDTmp*FactStep 
1863  
1864 
if (debug) WRITE(*,*) "DBG Main, L1475, FactStep=",FactStep 
1865 
if (debug.AND.(COORD.EQ.'ZMAT')) WRITE(*,'(A,12(1X,F10.6))') 'Scaled Step:',GradTmp(1:min(12,NFree)) 
1866  
1867 
IF (debug) THEN 
1868 
WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt 
1869 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1870 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1871 
ELSE 
1872 
DO Iat=1,Nat 
1873 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1874 
XyzGeomF(IGeom,1:3,Iat) 
1875 
END DO 
1876 
END IF 
1877 
END IF ! matches if (debug) THEN 
1878  
1879 
SELECT CASE (COORD) 
1880 
CASE ('ZMAT') 
1881 
IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) 
1882 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
1883 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
1884 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) 
1885 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
1886 
OrderInv(indzmat(2,2)),IntcoordF(IGeom,1) 
1887 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
1888 
OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi 
1889 
Idx=4 
1890 
DO Iat=4,Nat 
1891 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
1892 
OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & 
1893 
OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & 
1894 
OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi 
1895 
Idx=Idx+3 
1896 
END DO 
1897  
1898 
CASE ('BAKER') 
1899 
IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) 
1900 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
1901 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
1902 
WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull." 
1903  
1904 
CASE ('MIXED') 
1905 
IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) 
1906 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
1907 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
1908 
DO Iat=1,NCart 
1909 
WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat2:3*Iat) 
1910 
END DO 
1911 
Idx=3*NCart+1 
1912 
IF (NCart.EQ.1) THEN 
1913 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
1914 
OrderInv(indzmat(2,2)),IntcoordF(IGeom,4) 
1915 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
1916 
OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi 
1917  
1918 
Idx=7 
1919 
END IF 
1920 
IF (NCart.EQ.2) THEN 
1921 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
1922 
OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi 
1923 
Idx=9 
1924 
END IF 
1925  
1926 

1927 
DO Iat=max(NCart+1,4),Nat 
1928 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
1929 
OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & 
1930 
OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & 
1931 
OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi 
1932 
Idx=Idx+3 
1933 
END DO 
1934 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
1935  
1936 
CASE ('HYBRID','CART') 
1937 
XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
1938 
CASE DEFAULT 
1939 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP' 
1940 
STOP 
1941 
END SELECT 
1942  
1943 
IF (debug) THEN 
1944 
WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt 
1945 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1946 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1947 
ELSE 
1948 
DO Iat=1,Nat 
1949 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1950 
XyzGeomF(IGeom,1:3,Iat) 
1951 
END DO 
1952 
END IF 
1953 
END IF ! matches if (debug) THEN 
1954  
1955 
! We project out the tangential part of the gradient to know if we are done 
1956 
GradTmp=Grad(IGeom,:) 
1957 
NormGrad=sqrt(dot_product(GradTmp,GradTmp)) 
1958 
if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad 
1959 
SELECT CASE (COORD) 
1960 
CASE ('ZMAT','MIXED','BAKER') 
1961 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
1962 
Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:)) 
1963 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
1964 
Ca=S/(sqrt(Norm)*NormGrad) 
1965 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
1966 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
1967 
CASE ('CART','HYBRID') 
1968 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
1969 
Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:)) 
1970 
Ca=S/(sqrt(Norm)*NormGrad) 
1971 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
1972 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
1973 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
1974 
CASE DEFAULT 
1975 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP' 
1976 
STOP 
1977 
END SELECT 
1978 
IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi 
1979 
NormGrad=sqrt(DOT_Product(GradTmp,GradTmp)) 
1980 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
1981 
IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini 
1982  
1983 
END DO ! matches DO IGeom=2,NGeomF1 
1984 
! we save the old gradients 
1985 
GradOld=Grad 
1986 
EPathOld=Energies 
1987  
1988  
1989  
1990 
! Shall we reparameterize the path ? 
1991 
! PFL 2007/Apr/10 > 
1992 
! We call PathCreate in all cases, it will take care of the 
1993 
! reparameterization, as well as calculating the tangents. 
1994 
! IF (MOD(IOpt,IReparam).EQ.0) THEN 
1995 
XyzGeomI=XyzGeomF 
1996 
IntCoordI=IntCoordF 
1997  
1998 
Call PathCreate(IOpt) 
1999 
! END IF 
2000 
! < PFL 2007/Apr/10 
2001  
2002 
if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)') 'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini 
2003 
IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN 
2004  
2005 
! We have the new path, we calculate its energy and gradient 
2006  
2007 
Call EgradPath(IOpt) 
2008 
!IF(IOPT .EQ. 10) Then 
2009 
! Print *, 'Stopping at my checkpoint.' 
2010 
! STOP !This is my temporary checkpoint. 
2011 
!ENDIF 
2012  
2013 
! PFL 31 Mar 2008 > 
2014 
! Taken from v3.99 but we added a flag... 
2015 
! Added in v3.99 : MaxStep is modified according to the evolution of energies 
2016 
! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2017 
! else it is decreased by 0.8 
2018  
2019 
if (DynMaxStep) THEN 
2020 
if (debug) WRITE(*,*) "Dynamically updating the Maximum step" 
2021 
if (OptReac) THEN 
2022 
If (Energies(1)<EPathOld(1)) Then 
2023 
MaxStep(1)=min(1.1*MaxStep(1),2.*SMax) 
2024 
ELSE 
2025 
MaxStep(1)=max(0.8*MaxStep(1),SMax/2.) 
2026 
END IF 
2027 
END IF 
2028  
2029 
if (OptProd) THEN 
2030 
If (Energies(NGeomF)<EPathOld(NGeomF)) Then 
2031 
MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax) 
2032 
ELSE 
2033 
MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.) 
2034 
END IF 
2035 
END IF 
2036  
2037 
DO IGeom=2,NGeomF1 
2038 
If (Energies(IGeom)<EPathOld(IGeom)) Then 
2039 
MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax) 
2040 
ELSE 
2041 
MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.) 
2042 
END IF 
2043 
END DO 
2044 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
2045 
END IF ! test on DynMaxStep 
2046 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
2047 
! < PFL 31 MAr 2008 
2048 
! Also XyzGeomF is updated in EgradPath for Baker Case. 
2049 
! It should get updated for other cases too. 
2050  
2051 
DO IGeom=1,NGeomF 
2052 
SELECT CASE (COORD) 
2053 
CASE('ZMAT') 
2054 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2055 
GeomTmp=IntCoordF(IGeom,:) 
2056 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) 
2057 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
2058 
OrderInv(indzmat(2,2)),Geomtmp(1) 
2059 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
2060 
OrderInv(indzmat(3,2)),Geomtmp(2), & 
2061 
OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi 
2062 
Idx=4 
2063 
DO Iat=4,Nat 
2064 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
2065 
OrderInv(indzmat(iat,2)),Geomtmp(Idx), & 
2066 
OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, & 
2067 
OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi 
2068 
Idx=Idx+3 
2069 
END DO 
2070 
CASE('BAKER') 
2071 
!WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2072 
!GeomTmp=IntCoordF(IGeom,:) 
2073 
CASE('CART','HYBRID','MIXED') 
2074 
WRITE(IOOUT,'(1X,I5)') Nat 
2075 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2076 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2077 
DO I=1,Nat 
2078 
Iat=I 
2079 
If (renum) Iat=OrderInv(I) 
2080 
WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat2:3*Iat) 
2081 
END DO 
2082 
CASE DEFAULT 
2083 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP' 
2084 
STOP 
2085 
END SELECT 
2086 
END DO ! matches DO IGeom=1,NGeomF 
2087  
2088 
Call Write_path(IOpt) 
2089 
! Shall we analyze the geometries ? 
2090 
IF (AnaGeom) Call AnaPath(Iopt,IoDat) 
2091  
2092 
! We update the Hessian matrices 
2093 
IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord 
2094 
! First using the displacement from the old path 
2095 
IGeom0=2 
2096 
IGeomF=NGeomF1 
2097 
IF (OptReac) IGeom0=1 
2098 
IF (OptProd) IGeomF=NGeomF 
2099  
2100 
IF (FCalcHess) THEN 
2101 
DO IGeom=IGeom0,IGeomF 
2102 
SELECT CASE (COORD) 
2103 
CASE ('ZMAT','MIXED','BAKER') 
2104 
GeomTmp2=IntCoordF(IGeom,:) 
2105 
CASE ('CART','HYBRID') 
2106 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2107 
CASE DEFAULT 
2108 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
2109 
STOP 
2110 
END SELECT 
2111 
Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt) 
2112 
END DO 
2113 
ELSE 
2114 
DO IGeom=IGeom0,IGeomF 
2115 
GeomTmp=GeomOld_all(IGeom,:) 
2116 
SELECT CASE (COORD) 
2117 
CASE ('ZMAT','MIXED','BAKER') 
2118 
GeomTmp2=IntCoordF(IGeom,:) 
2119 
CASE ('CART','HYBRID') 
2120 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2121 
CASE DEFAULT 
2122 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
2123 
STOP 
2124 
END SELECT 
2125 
GeomTmp=GeomTmp2GeomTmp 
2126 
GradTmp=Grad(IGeom,:)GradOld(IGeom,:) 
2127 
HessTmp=Hess(IGeom,:,:) 
2128 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2129 
Hess(IGeom,:,:)=HessTmp 
2130 
END DO ! matches DO IGeom=IGeom0,IGeomF 
2131  
2132 
! Second using the neighbour points 
2133 
IF (HupNeighbour) THEN 
2134 
! Only one point for end points. 
2135 
IF (OptReac) THEN 
2136 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2137 
GeomTmp=IntCoordF(1,:)IntCoordF(2,:) 
2138 
ELSE 
2139 
GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))Reshape(XyzGeomF(2,:,:),(/3*Nat/)) 
2140 
END IF 
2141 
GradTmp=Grad(1,:)Grad(2,:) 
2142 
HessTmp=Hess(1,:,:) 
2143 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2144 
Hess(1,:,:)=HessTmp 
2145 
END IF 
2146  
2147 
IF (OptProd) THEN 
2148 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2149 
GeomTmp=IntCoordF(NGeomF,:)IntCoordF(NGeomF1,:) 
2150 
ELSE 
2151 
GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))Reshape(XyzGeomF(NGeomF1,:,:),(/3*Nat/)) 
2152 
END IF 
2153 
GradTmp=Grad(NGeomF,:)Grad(NGeomF1,:) 
2154 
HessTmp=Hess(NGeomF,:,:) 
2155 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2156 
Hess(NGeomF,:,:)=HessTmp 
2157 
END IF 
2158  
2159 
! Two points for the rest of the path 
2160 
DO IGeom=2,NGeomF1 
2161 
IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN 
2162 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom+1,:) 
2163 
ELSE 
2164 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) 
2165 
END IF 
2166 
GradTmp=Grad(IGeom,:)Grad(IGeom+1,:) 
2167 
HessTmp=Hess(IGeom,:,:) 
2168 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2169  
2170 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2171 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom1,:) 
2172 
ELSE 
2173 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom1,:,:),(/3*Nat/)) 
2174 
END IF 
2175 
GradTmp=Grad(IGeom,:)Grad(IGeom1,:) 
2176  
2177 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2178 
Hess(IGeom,:,:)=HessTmp 
2179 
END DO ! matches DO IGeom=2,NGeomF1 
2180 
END IF !matches IF (HupNeighbour) THEN 
2181 
END IF ! matches IF (FCalcHess) 
2182 
END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc) 
2183 
END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) 
2184  
2185 
! IF (PROG=="GAUSSIAN") THEN 
2186 
! DEALLOCATE(Gauss_paste) 
2187 
! DEALLOCATE(Gauss_root) 
2188 
! DEALLOCATE(Gauss_comment) 
2189 
! DEALLOCATE(Gauss_end) 
2190 
! END IF 
2191  
2192 
DEALLOCATE(XyzGeomF, IntCoordF) 
2193 
DEALLOCATE(XyzGeomI, IntCoordI) 
2194 
DEALLOCATE(XyzTangent,IntTangent) 
2195 
DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt) 
2196 
DEALLOCATE(GeomOld_all) 
2197  
2198 
if (PROG=="GAUSSIAN") THEN 
2199 
! We deallocate the Gauss_XX lists 
2200 
! transverse the list and deallocate each node 
2201 
previous => Gauss_root ! make current point to head of list 
2202 
DO 
2203 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2204 
current => previous%next ! make list point to next node of head 
2205 
DEALLOCATE(previous) ! deallocate current head node 
2206 
previous => current ! make current point to new head 
2207 
END DO 
2208  
2209 
previous => Gauss_comment ! make current point to head of list 
2210 
DO 
2211 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2212 
current => previous%next ! make list point to next node of head 
2213 
DEALLOCATE(previous) ! deallocate current head node 
2214 
previous => current ! make current point to new head 
2215 
END DO 
2216  
2217 
previous => Gauss_end ! make current point to head of list 
2218 
DO 
2219 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2220 
current => previous%next ! make list point to next node of head 
2221 
DEALLOCATE(previous) ! deallocate current head node 
2222 
previous => current ! make current point to new head 
2223 
END DO 
2224  
2225  
2226 
END IF 
2227  
2228 
DEALLOCATE(GeomTmp, Hess, GradTmp) 
2229 
IF (COORD.EQ.'ZMAT') DEALLOCATE(dzdc) 
2230 
IF (COORD.EQ.'BAKER') THEN 
2231 
DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t) 
2232 
DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local) 
2233 
END IF 
2234  
2235 
WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt 
2236  
2237 
Close(IOIN) 
2238 
Close(IOOUT) 
2239 
IF (AnaGeom) Close(IODAT) 
2240 
IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot) 
2241  
2242 
END PROGRAM PathOpt 