root / src / Path.f90 @ 9
History  View  Annotate  Download (89.5 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(*,*) "Opt'n 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 
! The following has been generated (March 19, 2013) using the Mini_help.tex file 
419 
! 1) Edit Mini_help.tex 
420 
! 2) latex Mini_help.tex 
421 
! 3) catdvi e 1 U Mini_help.dvi  sed re "s/\[U\+2022\]//g"  sed re "s/([^^[:space:]])\s+/\1 /g" > file.txt 
422 
! 4) edit file.txt to keep only the following lines, and to reformat a bit 
423 
! 5) awk '{print "WRITE(*,*) \"" $0 "\""}' file.txt > help 
424 
! 5) cut&paste help in Path.f90 ... 
425 
WRITE(*,*) " &path" 
426 
WRITE(*,*) " nat=3, ! Number of atoms" 
427 
WRITE(*,*) " ngeomi=3, ! Number of initial geometries" 
428 
WRITE(*,*) " ngeomf=12, !Number of geometries along the path" 
429 
WRITE(*,*) " OptReac=.T., ! Do you want to optimize the reactants ?" 
430 
WRITE(*,*) " OptProd=.T., ! Optimize the products" 
431 
WRITE(*,*) " coord='zmat', ! We use Zmatrix coordinates" 
432 
WRITE(*,*) " maxcyc=31, ! Max number of iterations" 
433 
WRITE(*,*) " IReparam=2,! redistribution of points along the path every 2 iterations" 
434 
WRITE(*,*) " ISpline=50, ! Start using spline interpolation at iteration 50" 
435 
WRITE(*,*) " Hinv=.T. , ! Use inverse of the Hessian internally (default: T)" 
436 
WRITE(*,*) " MW=T, ! Works in Mass Weighted coordiante (default T)" 
437 
WRITE(*,*) " PathName='Path_HCN_zmat_test', ! Name of the file used for path outputs" 
438 
WRITE(*,*) " prog='gaussian',! we use G03 to get energy and gradients" 
439 
WRITE(*,*) " SMax=0.1 ! Displacement cannot exceed 0.1 atomic units (or mass weighted at. unit)" 
440 
WRITE(*,*) " /" 
441 
WRITE(*,*) " 3" 
442 
WRITE(*,*) " Energy : 0.04937364" 
443 
WRITE(*,*) " H 0.0000 0.0000 0.0340" 
444 
WRITE(*,*) " C 0.0000 0.0000 1.1030" 
445 
WRITE(*,*) " N 0.0000 0.0000 2.2631" 
446 
WRITE(*,*) " 3" 
447 
WRITE(*,*) " Energy : 0.04937364" 
448 
WRITE(*,*) " H 0.0000 1.1000 1.1030" 
449 
WRITE(*,*) " C 0.0000 0.0000 1.1030" 
450 
WRITE(*,*) " N 0.0000 0.0000 2.2631" 
451 
WRITE(*,*) "3" 
452 
WRITE(*,*) " CNH" 
453 
WRITE(*,*) "H 0.000000 0.000000 3.3" 
454 
WRITE(*,*) "C 0.000000 0.000000 1.1" 
455 
WRITE(*,*) "N 0.000000 0.000000 2.26" 
456 
WRITE(*,*) "%chk=Test" 
457 
WRITE(*,*) "#P AM1 FORCE" 
458 
WRITE(*,*) " HCN est bien" 
459 
WRITE(*,*) "" 
460 
WRITE(*,*) "0,1" 
461 
WRITE(*,*) "H 0.000000 0.000000 0.000000" 
462 
WRITE(*,*) "C 0.000000 0.000000 1.000" 
463 
WRITE(*,*) "N 0.000000 0.000000 3.00" 
464 
WRITE(*,*) "" 
465 
WRITE(*,*) "*******" 
466 
WRITE(*,*) "* Compulsory variables are:" 
467 
WRITE(*,*) "*******" 
468 
WRITE(*,*) "NGeomi: Number of geometries defining the Initial path" 
469 
WRITE(*,*) "NGeomf: Number of geometries defining the Final path" 
470 
WRITE(*,*) "Nat: Number of atoms" 
471 
WRITE(*,*) "" 
472 
WRITE(*,*) "*******" 
473 
WRITE(*,*) "* Other options:" 
474 
WRITE(*,*) "*******" 
475 
WRITE(*,*) "Input: String that indicates the type of the input geometries. Accepted values are: Cart (or" 
476 
WRITE(*,*) " Xmol or Xyz), Vasp, Turbomole or Siesta." 
477 
WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations." 
478 
WRITE(*,*) " Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext." 
479 
WRITE(*,*) "" 
480 
WRITE(*,*) "" 
481 
WRITE(*,*) "  In case of a Gaussian calculations, input must be set to Cart. One example of a gaus" 
482 
WRITE(*,*) " sian input should be added at the end of the input file.See example file Test_HCN_zmat_g03.path." 
483 
WRITE(*,*) "  In the case of a VASP calculation, if input is set to Cart, then the preamble of a" 
484 
WRITE(*,*) " VASP calculation s" 
485 
WRITE(*,*) "  Mopac: Examples using sequential call to MOPAC2009 to calculate the energies and" 
486 
WRITE(*,*) " forces along the path. hould be added at the end of the input file. See example file" 
487 
WRITE(*,*) " Test_VASP_cart.path. In the case of a VASP calculation, one should also give value" 
488 
WRITE(*,*) " of the RunMode variable ." 
489 
WRITE(*,*) "  In the case of a SIESTA calculation, an example of a Siesta input file should be added" 
490 
WRITE(*,*) " at the end of the input file. See Test_Siesta.path." 
491 
WRITE(*,*) "" 
492 
WRITE(*,*) "RunMode: When running on a multiprocessor machine, this indicates wether Opt'n Path" 
493 
WRITE(*,*) " should calculate the energy and gradient of the whole path in parallel or not. User has" 
494 
WRITE(*,*) " two options. One is to calculate the energy and gradient of each point sequentially. This" 
495 
WRITE(*,*) " is usefull when running on one (or two) processors. In this case, RunMode should be put" 
496 
WRITE(*,*) " to SERIAL.When running in parallel with 8 or more processors, one can use VASP to" 
497 
WRITE(*,*) " calculate simultaneously the energies and gradients for all points, as in a normal NEB cal" 
498 
WRITE(*,*) " culation. In this case, RunMode must be set to PARA. For now, this is usefull only for VASP." 
499 
WRITE(*,*) "" 
500 
WRITE(*,*) "ProgExe: Name (with full path) of the executable to be used to get energies and gradients." 
501 
WRITE(*,*) " For example, if VASP is used in parallel, one might have something like:" 
502 
WRITE(*,*) " ProgExe='/usr/local/mpich/bin/mpirun machinefile machine np 8 ~/bin/VASP_46'." 
503 
WRITE(*,*) " Another option that I use, is to put ProgExe='./run_vasp' and I put every option needed" 
504 
WRITE(*,*) " to run VASP into the run_vasp file." 
505 
WRITE(*,*) "EReac: (REAL) By default, Opt'n Path does not compute the energy of the reactants and" 
506 
WRITE(*,*) " products. This thus indicates the reactants energy to Opt'n Path to have better plots at" 
507 
WRITE(*,*) " the end." 
508 
WRITE(*,*) "EProd: (REAL) By default, Opt'n Path does not compute the energy of the reactants and" 
509 
WRITE(*,*) " products. This thus indicates the products energy to Opt'n Path to have better plots." 
510 
WRITE(*,*) "" 
511 
WRITE(*,*) "PathName: Prefix used to save the path. Default is Path." 
512 
WRITE(*,*) "Poscar: string that will be used as the prefix for the different POSCAR files in a VASP calcu" 
513 
WRITE(*,*) " lations. Usefull only if PathOnly=.TRUE., not used for internal calculations." 
514 
WRITE(*,*) "CalcName: Prefix for the files used for the energy and gradient calculations" 
515 
WRITE(*,*) "" 
516 
WRITE(*,*) "ISuffix: Suffix for the input file used for energy and gradient calculations. The full inputfile" 
517 
WRITE(*,*) " name will thus be CalcName.ISuffix." 
518 
WRITE(*,*) "OSuffix: Suffix for the output file used for energy and gradient calculations. The full outputfile" 
519 
WRITE(*,*) " name will thus be CalcName.OSuffix." 
520 
WRITE(*,*) "" 
521 
WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates. Valid only for" 
522 
WRITE(*,*) " Coord=Zmat, Hybrid or Mixed." 
523 
WRITE(*,*) "Fact: REAL used to define if two atoms are linked. If d(A,B) =< fact* (rcov(A) + rcov(B))," 
524 
WRITE(*,*) " then A and B are considered Linked." 
525 
WRITE(*,*) "" 
526 
WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info." 
527 
WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:" 
528 
WRITE(*,*) "  CART (or Xyz): works in cartesian" 
529 
WRITE(*,*) "  Zmat: works in internal coordinates (Zmat)" 
530 
WRITE(*,*) "  Mixed: frozen atoms, as well as atoms defined by the 'cart' array(see below) are" 
531 
WRITE(*,*) " describe in CARTESIAN, whereas the others are described in Zmat" 
532 
WRITE(*,*) "  Baker: use of Baker coordinates, also called delocalized internal coordinates" 
533 
WRITE(*,*) "  Hybrid: geometries are described in zmat, but the gradient are used in cartesian" 
534 
WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
535 
WRITE(*,*) "" 
536 
WRITE(*,*) "  RFO: Rational function optimization" 
537 
WRITE(*,*) "  GDIIS: Geometry optimization using direct inversion in the iterative supspace" 
538 
WRITE(*,*) "HUpdate: method to update the hessian. By default, it is MurtaghSargent Except for geom" 
539 
WRITE(*,*) " etry optimization where it is BFGS." 
540 
WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization" 
541 
WRITE(*,*) "" 
542 
WRITE(*,*) "Smax: Maximum length of a step during path optimization" 
543 
WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary" 
544 
WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is" 
545 
WRITE(*,*) " taken" 
546 
WRITE(*,*) "" 
547 
WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path" 
548 
WRITE(*,*) " that is kept. FTan=1. corresponds to the full displacement, whereas FTan=0. gives a" 
549 
WRITE(*,*) " displacement orthogonal to the path." 
550 
WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of" 
551 
WRITE(*,*) " reparameterization." 
552 
WRITE(*,*) "IReparamT: When the path is not reparameterised at each iteration, this gives the frequency" 
553 
WRITE(*,*) " of reparameterization of the tangents." 
554 
WRITE(*,*) "" 
555 
WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path. This option indicates" 
556 
WRITE(*,*) " the first step where spline interpolation is used." 
557 
WRITE(*,*) "BoxTol: Real between 0. and 1. When doing periodic calculations, it might happen that an" 
558 
WRITE(*,*) " atom moves out of the unit cell. Path detects this by comparing the displacement to" 
559 
WRITE(*,*) " boxtol: if an atom moves by more than Boxtol, then it is moved back into the unit cell." 
560 
WRITE(*,*) " Default value: 0.5." 
561 
WRITE(*,*) "FrozTol: (Real) This indicates the threshold to determine wether an atom moves between two" 
562 
WRITE(*,*) " images. Default is 1e4." 
563 
WRITE(*,*) "" 
564 
WRITE(*,*) "OptGeom: This INTEGER indicates the index of the geometry you want to optimize. If" 
565 
WRITE(*,*) " OptGeom is set, then Opt'n Path performs a geometry optimization instead of a path" 
566 
WRITE(*,*) " interpolation." 
567 
WRITE(*,*) "GeomFile: Name of the file to print the geometries and their energy during a geometry opti" 
568 
WRITE(*,*) " mization. If this variable is not given then nothing is printed." 
569 
WRITE(*,*) "AnaFile: Name of the file to print the values of the geometrical parameters that are monitored" 
570 
WRITE(*,*) " if AnaGeom=.TRUE.. Default is PathName.dat" 
571 
WRITE(*,*) "" 
572 
WRITE(*,*) "GplotFile: Name of the gnuplot file to plot the evolution of the geometrical parameters that" 
573 
WRITE(*,*) " are monitored if AnaGeom=.TRUE.. Default is PathName.gplot" 
574 
WRITE(*,*) "" 
575 
WRITE(*,*) "*******" 
576 
WRITE(*,*) "* Arrays:" 
577 
WRITE(*,*) "*******" 
578 
WRITE(*,*) "" 
579 
WRITE(*,*) "Rcov: Array containing the covalent radii of the first 86 elements. You can modify it using," 
580 
WRITE(*,*) " rcov(6)=0.8." 
581 
WRITE(*,*) "Mass: Array containing the atomic mass of the first 86 elements." 
582 
WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations. If not given, Opt'n Path will" 
583 
WRITE(*,*) " read the POTCAR file." 
584 
WRITE(*,*) "" 
585 
WRITE(*,*) "*******" 
586 
WRITE(*,*) "* Flags:" 
587 
WRITE(*,*) "*******" 
588 
WRITE(*,*) "" 
589 
WRITE(*,*) "MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE." 
590 
WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE." 
591 
WRITE(*,*) " unless for Coord equals CART." 
592 
WRITE(*,*) "" 
593 
WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products." 
594 
WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants." 
595 
WRITE(*,*) "CalcEProd: if TRUE the product energy will be computed. Default is FALSE. Not Compatible" 
596 
WRITE(*,*) " with RunMode=Para." 
597 
WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compat" 
598 
WRITE(*,*) " ible with RunMode=Para." 
599 
WRITE(*,*) "" 
600 
WRITE(*,*) "PathOnly: TRUE if one wants to generate the initial path, and stops." 
601 
WRITE(*,*) "Align: If .FALSE., successive geometries along the path are not aligned on each other before" 
602 
WRITE(*,*) " path interpolation. Default is .TRUE. if there are 4 atoms or more." 
603 
WRITE(*,*) "Hinv: if True, then Hessian inversed is used." 
604 
WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations." 
605 
WRITE(*,*) "" 
606 
WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of" 
607 
WRITE(*,*) " the path." 
608 
WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms. If True, a &frozenlist" 
609 
WRITE(*,*) " namelist containing the list of frozen atoms must be given. If VASP is used, and frozen is" 
610 
WRITE(*,*) " not given here, the program will use the F flags of the input geometry" 
611 
WRITE(*,*) "FCart: True if one wants to describe some atoms using cartesian coordinates. *** Only used in" 
612 
WRITE(*,*) " 'mixed' calculations. *** If True, a &cartlist namelist containing the list of cart atoms" 
613 
WRITE(*,*) " must be given. By default, only frozen atoms are described in cartesian coordinates." 
614 
WRITE(*,*) "" 
615 
WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms." 
616 
WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now." 
617 
WRITE(*,*) "" 
618 
WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files. See also" 
619 
WRITE(*,*) " the POSCAR option. This can be used only if prog or input=VASP." 
620 
WRITE(*,*) "AnaGeom: If TRUE, Opt'n Path will create a file .dat for geometries analysis. If True, " 
621 
WRITE(*,*) "Opt'n Path will look for the AnaList namelist after the Path Namelist." 
622 
WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry." 
623 
WRITE(*,*) " If energy goes up, Smax = Smax*0.8, if not Smax = Smax*1.2. It is ensured that the" 
624 
WRITE(*,*) " dynamical Smax is within [0.5*SMax_0 ,2*Smax_0 ]." 
625 
WRITE(*,*) "" 
626 
WRITE(*,*) "*******" 
627 
WRITE(*,*) "* Additional namelists:" 
628 
WRITE(*,*) "*******" 
629 
WRITE(*,*) "" 
630 
WRITE(*,*) "&cartlist list= ... &end This gives the list of atoms that are described using cartesian coor" 
631 
WRITE(*,*) " dinates. Read only if FCart=.TRUE.. To indicate an atom range, from 1 to 5 for example," 
632 
WRITE(*,*) " one can put 1 5 in the list. For example:" 
633 
WRITE(*,*) " &cartlist list = 1 2 6 12 20 &end" 
634 
WRITE(*,*) " will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian." 
635 
WRITE(*,*) "&Frozenlist list= ... &end This gives the list of atoms that are frozen during optimization." 
636 
WRITE(*,*) " Read only if FFrozen=.TRUE.. To indicate an atom range, from 1 to 5 for example, one" 
637 
WRITE(*,*) " can put 1 5 in the list." 
638 
WRITE(*,*) "" 
639 
WRITE(*,*) "&Analist nb= ... &end list of variables for geometry analysis. If present and if Ana" 
640 
WRITE(*,*) " Geom=T then Opt'n Path will read it and print the values of the variable in a .dat file If" 
641 
WRITE(*,*) " AnaGeom is T but Analist is not given, then Path.dat just contains the energy. Analist" 
642 
WRITE(*,*) " contains the number of variables to monitor: Nb, and is followed by the description of the" 
643 
WRITE(*,*) " variables among: b(ond) At1 At2 a(ngle) At1 At2 At3 d(ihedral) At1 At2 At3 At4 c NbAt" 
644 
WRITE(*,*) " At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added at the" 
645 
WRITE(*,*) " end of the real atoms of the system." 
646 
!!!!!!!!!! end of included file 
647 
WRITE(*,*) "" 
648 
STOP 
649 
END IF 
650  
651 
OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') 
652 
IOOUT=6 
653 
CASE (0) 
654 
IOIN=5 
655 
IOOUT=6 
656 
END SELECT 
657  
658  
659 
! We initiliaze variables 
660 
Pi=dacos(1.0d0) 
661 
Nat=1 
662 
Fact=1.1 
663 
IGeomRef=1 
664 
NGeomI=1 
665 
NGeomF=8 
666 
IReparam=5 
667 
IReparamT=1 
668 
ISpline=5 
669 
debug=.False. 
670 
MW=.TRUE. 
671 
bohr=.False. 
672 
Coord='ZMAT' 
673 
Step_method='RFO' 
674 
DebugFile='Path.valid' 
675 
PathName="Path" 
676 
MaxCyc=50 
677 
SMax=0.1D0 
678 
SThresh=1. 
679 
GThresh=1. 
680 
FTan=0. 
681 
Hinv=.TRUE. 
682 
IniHUp=.FALSE. 
683 
HupNeighbour=.FALSE. 
684 
HesUpd="EMPTY" 
685 
HUpdate="EMPTY" 
686 
FFrozen=.FALSE. 
687 
FCart=.FALSE. 
688 
List=0 
689 
renum=.TRUE. 
690 
OptReac=.FALSE. 
691 
OptProd=.FALSE. 
692 
CalcEReac=.FALSE. 
693 
CalcEProd=.FALSE. 
694 
EReac=0.d0 
695 
EProd=0.d0 
696 
OptGeom=1 
697 
GeomFile="EMPTY" 
698 
AnaGeom=.FALSE. 
699 
AnaFile="EMPTY" 
700 
GplotFile="EMPTY" 
701 
Nb=0 
702 
NbCom=0 
703 
PathOnly=.FALSE. 
704 
Prog='EMPTY' 
705 
ProgExe='EMPTY' 
706 
Runmode='SERIAL' 
707 
Input='EMPTY' 
708 
Poscar="POSCAR" 
709 
AutoCart=.FALSE. 
710 
VMD=.TRUE. 
711 
WriteVASP=.FALSE. 
712 
DynMaxStep=.FALSE. 
713 
CalcName="EMPTY" 
714 
ISuffix="EMPTY" 
715 
OSuffix="EMPTY" 
716 
NGintMax=10 
717 
Align=.TRUE. 
718  
719 
! Number of point used to compute the length of the path, 
720 
! and to reparameterize the path 
721 
NMaxPtPath=1 
722 
FrozTol=1e4 
723  
724 
! Read the namelist 
725 
read (IOIN,path) 
726  
727 
debug=valid("Main").or.valid("Path") 
728  
729 
FCalcHess=valid("CalcHess") 
730 
PathName=AdjustL(PathName) 
731 
Coord=AdjustL(Coord) 
732 
CALL UpCase(coord) 
733  
734 
Step_method=AdjustL(Step_method) 
735 
CALL UpCase(Step_method) 
736  
737 
Prog=AdjustL(Prog) 
738 
CALL UpCase(prog) 
739  
740 
Input=AdjustL(Input) 
741 
CALL UpCase(input) 
742  
743 
Runmode=AdjustL(Runmode) 
744 
Call UpCase(Runmode) 
745 
IF (Runmode(1:4).EQ.'PARA') Runmode="PARA" 
746  
747 
if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN 
748 
WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'" 
749 
WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE" 
750 
CalcEReac=.FALSE. 
751 
CalcEProd=.FALSE. 
752 
END IF 
753  
754 
If (IReparamT<=0) IReparamT=IReparam 
755  
756 
! We take care of HesUpd flag 
757 
! this is needed because some users might still use it 
758 
if (HesUpd.NE."EMPTY") THEN 
759 
! We are pityless: 
760 
WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead" 
761 
WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead" 
762 
STOP 
763 
END IF 
764  
765 
IF (HUpdate=='EMPTY') THEN 
766 
IF (OptGeom>=1) THEN 
767 
HupDate="BFGS" 
768 
ELSE 
769 
HUpdate="MS" 
770 
END IF 
771 
END IF 
772  
773 
If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 
774 
COORD='CART' 
775 
END IF 
776  
777 
IF (COORD.EQ.'CART') THEN 
778 
Renum=.FALSE. 
779 
IF (OptGeom>0) THEN 
780 
IGeomRef=OptGeom 
781 
ELSE 
782 
IGeomRef=1 
783 
END IF 
784 
END IF 
785  
786 

787 
if (Prog.EQ.'TEST2') Prog="CHAMFRE" 
788 
if (Prog.EQ.'2D') Prog="TEST2D" 
789 
if (Prog.EQ.'TEST3') Prog="LEPS" 
790  
791  
792 
Select Case (Prog) 
793 
CASE ("EMPTY","G09","GAUSSIAN") 
794 
Prog='GAUSSIAN' 
795 
If (ProgExe=="EMPTY") ProgExe="g09" 
796 
if (CalcName=="EMPTY") CalcName="Path" 
797 
if (ISuffix=="EMPTY") ISuffix="com" 
798 
if (OSuffix=="EMPTY") OSuffix="log" 
799 
UnitProg="au" 
800 
ConvE=au2kcal 
801 
CASE ("G03") 
802 
Prog='GAUSSIAN' 
803 
If (ProgExe=="EMPTY") ProgExe="g03" 
804 
if (CalcName=="EMPTY") CalcName="Path" 
805 
if (ISuffix=="EMPTY") ISuffix="com" 
806 
if (OSuffix=="EMPTY") OSuffix="log" 
807 
UnitProg="au" 
808 
ConvE=au2kcal 
809 
CASE ("EXT") 
810 
If (ProgExe=="EMPTY") ProgExe="./Ext.exe" 
811 
if (CalcName=="EMPTY") CalcName="Ext" 
812 
if (ISuffix=="EMPTY") ISuffix="in" 
813 
if (OSuffix=="EMPTY") OSuffix="out" 
814 
UnitProg="au" 
815 
ConvE=au2kcal 
816 
CASE ('MOPAC') 
817 
If (ProgExe=="EMPTY") ProgExe="mopac" 
818 
if (CalcName=="EMPTY") CalcName="Path" 
819 
if (ISuffix=="EMPTY") ISuffix="mop" 
820 
if (OSuffix=="EMPTY") OSuffix="out" 
821 
UnitProg="au" 
822 
ConvE=au2kcal 
823 
CASE ("SIESTA") 
824 
If (ProgExe=="EMPTY") ProgExe="siesta" 
825 
if (CalcName=="EMPTY") CalcName="Path" 
826 
if (ISuffix=="EMPTY") ISuffix="fdf" 
827 
if (OSuffix=="EMPTY") OSuffix="out" 
828 
UnitProg="eV" 
829 
ConvE=eV2kcal 
830 
CASE ("VASP") 
831 
If (ProgExe=="EMPTY") ProgExe="VASP" 
832 
UnitProg="eV" 
833 
ConvE=eV2kcal 
834 
CASE ("TM","TURBOMOLE") 
835 
Prog="TURBOMOLE" 
836 
If (ProgExe=="EMPTY") ProgExe="jobex c 1 keep" 
837 
UnitProg="au" 
838 
ConvE=au2kcal 
839 
CASE ("TEST","CHAMFRE","TEST2D","LEPS") 
840 
ProgExe="" 
841 
UnitProg="au" 
842 
ConvE=au2kcal 
843 
CASE DEFAULT 
844 
WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP" 
845 
STOP 
846 
END SELECT 
847  
848 
if (Input.EQ.'EMPTY') THEN 
849 
Select Case (Prog) 
850 
CASE ("VASP") 
851 
Input=Prog 
852 
CASE ("CHAMFRE") 
853 
Input="VASP" 
854 
CASE DEFAULT 
855 
Input='XYZ' 
856 
END SELECT 
857 
WRITE(*,*) "Input has been set to the default: ",INPUT 
858 
END IF 
859  
860 
WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) 
861  
862 
IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN 
863 
WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode." 
864 
STOP 
865 
END IF 
866  
867 
if (OptGeom.GE.1) THEN 
868 
FPrintGeom=.FALSE. 
869 
If (GeomFile/='EMPTY') THEN 
870 
FPrintGeom=.TRUE. 
871 
END IF 
872 
IF (IGeomRef.LE.0) THEN 
873 
IGeomRef=OptGeom 
874 
ELSE 
875 
IF (IGeomRef.NE.OptGeom) THEN 
876 
WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
877 
WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
878  
879 
WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
880 
WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
881 
STOP 
882 
END IF 
883 
END IF 
884 
END IF 
885  
886 
ALLOCATE(Cart(Nat)) 
887 
Cart=0 
888  
889 
ALLOCATE(Frozen(Nat)) 
890 
Frozen=0 
891  
892 
IF (FCart) THEN 
893 
! We rewind the file to be sure to browse all of it to read the Cart namelist 
894 
REWIND(IOIN) 
895 
List=0 
896 
READ(IOIN,cartlist) 
897 
IF (COORD.NE.'CART') THEN 
898 
Cart=List(1:Nat) 
899 
if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart 
900 
ELSE 
901 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
902 
WRITE(*,*) "I will discard the cartlist" 
903 
Cart=0 
904 
END IF 
905 
END IF 
906  
907 
if (FFrozen) THEN 
908  
909 
REWIND(IOIN) 
910 
List=0 
911 
READ(IOIN,Frozenlist) 
912 
Frozen=List(1:Nat) 
913 
if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen 
914 
END IF 
915  
916 
If (FCart) Call Expand(Cart,NCart,Nat) 
917 
IF (FFrozen) Call Expand(Frozen,NFroz,Nat) 
918  
919 
if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart 
920 
if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen 
921  
922 
if (AnaGeom) THEN 
923 
REWIND(IOIN) 
924 
READ(IOIN,AnaList,IOSTAT=IoS) 
925 
IF (IOS==0) THEN ! we have a AnaList namelist 
926 
Call ReadAnaList 
927 
ELSE 
928 
! No AnaList namelist, we will print only the energy 
929 
FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)' 
930 
END IF 
931 
IF (AnaFile=="EMPTY") THEN 
932 
AnaFile=TRIM(PathName) // '.dat' 
933 
END IF 
934 
OPEN(IODat,File=AnaFile) 
935 
IF (GplotFile=="EMPTY") THEN 
936 
GplotFile=TRIM(PathName) // '.gplot' 
937 
END IF 
938 

939 
OPEN(IODat,File=AnaFile) 
940 
If (OptGeom<=0) THEN 
941 
OPEN(IOGplot,File=GPlotFile) 
942 
WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot persist" 
943 
WRITE(IoGplot,'(A)') " set pointsize 2" 
944 
WRITE(IoGplot,'(A)') " xcol=1" 
945 
WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3 
946 
END IF 
947 
END IF 
948  
949 

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

2040 
DO Iat=max(NCart+1,4),Nat 
2041 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
2042 
OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & 
2043 
OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & 
2044 
OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi 
2045 
Idx=Idx+3 
2046 
END DO 
2047 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
2048  
2049 
CASE ('HYBRID','CART') 
2050 
XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
2051 
CASE DEFAULT 
2052 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP' 
2053 
STOP 
2054 
END SELECT 
2055  
2056 
IF (debug) THEN 
2057 
WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt 
2058 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2059 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
2060 
ELSE 
2061 
DO Iat=1,Nat 
2062 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
2063 
XyzGeomF(IGeom,1:3,Iat) 
2064 
END DO 
2065 
END IF 
2066 
END IF ! matches if (debug) THEN 
2067  
2068 
! We project out the tangential part of the gradient to know if we are done 
2069 
GradTmp=Grad(IGeom,:) 
2070 
NormGrad=sqrt(dot_product(GradTmp,GradTmp)) 
2071 
if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad 
2072 
SELECT CASE (COORD) 
2073 
CASE ('ZMAT','MIXED','BAKER') 
2074 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
2075 
Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:)) 
2076 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
2077 
Ca=S/(sqrt(Norm)*NormGrad) 
2078 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
2079 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
2080 
CASE ('CART','HYBRID') 
2081 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
2082 
Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:)) 
2083 
Ca=S/(sqrt(Norm)*NormGrad) 
2084 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
2085 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
2086 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
2087 
CASE DEFAULT 
2088 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP' 
2089 
STOP 
2090 
END SELECT 
2091 
IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi 
2092 
NormGrad=sqrt(DOT_Product(GradTmp,GradTmp)) 
2093 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
2094 
IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini 
2095  
2096 
END DO ! matches DO IGeom=2,NGeomF1 
2097 
! we save the old gradients 
2098 
GradOld=Grad 
2099 
EPathOld=Energies 
2100  
2101  
2102  
2103 
! Shall we reparameterize the path ? 
2104 
! PFL 2007/Apr/10 > 
2105 
! We call PathCreate in all cases, it will take care of the 
2106 
! reparameterization, as well as calculating the tangents. 
2107 
! IF (MOD(IOpt,IReparam).EQ.0) THEN 
2108 
XyzGeomI=XyzGeomF 
2109 
IntCoordI=IntCoordF 
2110  
2111 
Call PathCreate(IOpt) 
2112 
! END IF 
2113 
! < PFL 2007/Apr/10 
2114  
2115 
if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)') 'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini 
2116 
IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN 
2117  
2118 
! We have the new path, we calculate its energy and gradient 
2119  
2120 
Call EgradPath(IOpt) 
2121 
!IF(IOPT .EQ. 10) Then 
2122 
! Print *, 'Stopping at my checkpoint.' 
2123 
! STOP !This is my temporary checkpoint. 
2124 
!ENDIF 
2125  
2126 
! PFL 31 Mar 2008 > 
2127 
! Taken from v3.99 but we added a flag... 
2128 
! Added in v3.99 : MaxStep is modified according to the evolution of energies 
2129 
! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2130 
! else it is decreased by 0.8 
2131  
2132 
if (DynMaxStep) THEN 
2133 
if (debug) WRITE(*,*) "Dynamically updating the Maximum step" 
2134 
if (OptReac) THEN 
2135 
If (Energies(1)<EPathOld(1)) Then 
2136 
MaxStep(1)=min(1.1*MaxStep(1),2.*SMax) 
2137 
ELSE 
2138 
MaxStep(1)=max(0.8*MaxStep(1),SMax/2.) 
2139 
END IF 
2140 
END IF 
2141  
2142 
if (OptProd) THEN 
2143 
If (Energies(NGeomF)<EPathOld(NGeomF)) Then 
2144 
MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax) 
2145 
ELSE 
2146 
MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.) 
2147 
END IF 
2148 
END IF 
2149  
2150 
DO IGeom=2,NGeomF1 
2151 
If (Energies(IGeom)<EPathOld(IGeom)) Then 
2152 
MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax) 
2153 
ELSE 
2154 
MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.) 
2155 
END IF 
2156 
END DO 
2157 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
2158 
END IF ! test on DynMaxStep 
2159 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
2160 
! < PFL 31 MAr 2008 
2161 
! Also XyzGeomF is updated in EgradPath for Baker Case. 
2162 
! It should get updated for other cases too. 
2163  
2164 
DO IGeom=1,NGeomF 
2165 
SELECT CASE (COORD) 
2166 
CASE('ZMAT') 
2167 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2168 
GeomTmp=IntCoordF(IGeom,:) 
2169 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) 
2170 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
2171 
OrderInv(indzmat(2,2)),Geomtmp(1) 
2172 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
2173 
OrderInv(indzmat(3,2)),Geomtmp(2), & 
2174 
OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi 
2175 
Idx=4 
2176 
DO Iat=4,Nat 
2177 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
2178 
OrderInv(indzmat(iat,2)),Geomtmp(Idx), & 
2179 
OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, & 
2180 
OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi 
2181 
Idx=Idx+3 
2182 
END DO 
2183 
CASE('BAKER') 
2184 
!WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2185 
!GeomTmp=IntCoordF(IGeom,:) 
2186 
CASE('CART','HYBRID','MIXED') 
2187 
WRITE(IOOUT,'(1X,I5)') Nat 
2188 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2189 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2190 
DO I=1,Nat 
2191 
Iat=I 
2192 
If (renum) Iat=OrderInv(I) 
2193 
WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat2:3*Iat) 
2194 
END DO 
2195 
CASE DEFAULT 
2196 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP' 
2197 
STOP 
2198 
END SELECT 
2199 
END DO ! matches DO IGeom=1,NGeomF 
2200  
2201 
Call Write_path(IOpt) 
2202 
! Shall we analyze the geometries ? 
2203 
IF (AnaGeom) Call AnaPath(Iopt,IoDat) 
2204  
2205 
! We update the Hessian matrices 
2206 
IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord 
2207 
! First using the displacement from the old path 
2208 
IGeom0=2 
2209 
IGeomF=NGeomF1 
2210 
IF (OptReac) IGeom0=1 
2211 
IF (OptProd) IGeomF=NGeomF 
2212  
2213 
IF (FCalcHess) THEN 
2214 
DO IGeom=IGeom0,IGeomF 
2215 
SELECT CASE (COORD) 
2216 
CASE ('ZMAT','MIXED','BAKER') 
2217 
GeomTmp2=IntCoordF(IGeom,:) 
2218 
CASE ('CART','HYBRID') 
2219 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2220 
CASE DEFAULT 
2221 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
2222 
STOP 
2223 
END SELECT 
2224 
Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt) 
2225 
END DO 
2226 
ELSE 
2227 
DO IGeom=IGeom0,IGeomF 
2228 
GeomTmp=GeomOld_all(IGeom,:) 
2229 
SELECT CASE (COORD) 
2230 
CASE ('ZMAT','MIXED','BAKER') 
2231 
GeomTmp2=IntCoordF(IGeom,:) 
2232 
CASE ('CART','HYBRID') 
2233 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2234 
CASE DEFAULT 
2235 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
2236 
STOP 
2237 
END SELECT 
2238 
GeomTmp=GeomTmp2GeomTmp 
2239 
GradTmp=Grad(IGeom,:)GradOld(IGeom,:) 
2240 
HessTmp=Hess(IGeom,:,:) 
2241 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2242 
Hess(IGeom,:,:)=HessTmp 
2243 
END DO ! matches DO IGeom=IGeom0,IGeomF 
2244  
2245 
! Second using the neighbour points 
2246 
IF (HupNeighbour) THEN 
2247 
! Only one point for end points. 
2248 
IF (OptReac) THEN 
2249 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2250 
GeomTmp=IntCoordF(1,:)IntCoordF(2,:) 
2251 
ELSE 
2252 
GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))Reshape(XyzGeomF(2,:,:),(/3*Nat/)) 
2253 
END IF 
2254 
GradTmp=Grad(1,:)Grad(2,:) 
2255 
HessTmp=Hess(1,:,:) 
2256 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2257 
Hess(1,:,:)=HessTmp 
2258 
END IF 
2259  
2260 
IF (OptProd) THEN 
2261 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2262 
GeomTmp=IntCoordF(NGeomF,:)IntCoordF(NGeomF1,:) 
2263 
ELSE 
2264 
GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))Reshape(XyzGeomF(NGeomF1,:,:),(/3*Nat/)) 
2265 
END IF 
2266 
GradTmp=Grad(NGeomF,:)Grad(NGeomF1,:) 
2267 
HessTmp=Hess(NGeomF,:,:) 
2268 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2269 
Hess(NGeomF,:,:)=HessTmp 
2270 
END IF 
2271  
2272 
! Two points for the rest of the path 
2273 
DO IGeom=2,NGeomF1 
2274 
IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN 
2275 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom+1,:) 
2276 
ELSE 
2277 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) 
2278 
END IF 
2279 
GradTmp=Grad(IGeom,:)Grad(IGeom+1,:) 
2280 
HessTmp=Hess(IGeom,:,:) 
2281 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2282  
2283 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2284 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom1,:) 
2285 
ELSE 
2286 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom1,:,:),(/3*Nat/)) 
2287 
END IF 
2288 
GradTmp=Grad(IGeom,:)Grad(IGeom1,:) 
2289  
2290 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2291 
Hess(IGeom,:,:)=HessTmp 
2292 
END DO ! matches DO IGeom=2,NGeomF1 
2293 
END IF !matches IF (HupNeighbour) THEN 
2294 
END IF ! matches IF (FCalcHess) 
2295 
END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc) 
2296 
END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) 
2297  
2298 
! IF (PROG=="GAUSSIAN") THEN 
2299 
! DEALLOCATE(Gauss_paste) 
2300 
! DEALLOCATE(Gauss_root) 
2301 
! DEALLOCATE(Gauss_comment) 
2302 
! DEALLOCATE(Gauss_end) 
2303 
! END IF 
2304  
2305 
DEALLOCATE(XyzGeomF, IntCoordF) 
2306 
DEALLOCATE(XyzGeomI, IntCoordI) 
2307 
DEALLOCATE(XyzTangent,IntTangent) 
2308 
DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt) 
2309 
DEALLOCATE(GeomOld_all) 
2310  
2311 
if (PROG=="GAUSSIAN") THEN 
2312 
! We deallocate the Gauss_XX lists 
2313 
! transverse the list and deallocate each node 
2314 
previous => Gauss_root ! make current point to head of list 
2315 
DO 
2316 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2317 
current => previous%next ! make list point to next node of head 
2318 
DEALLOCATE(previous) ! deallocate current head node 
2319 
previous => current ! make current point to new head 
2320 
END DO 
2321  
2322 
previous => Gauss_comment ! make current point to head of list 
2323 
DO 
2324 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2325 
current => previous%next ! make list point to next node of head 
2326 
DEALLOCATE(previous) ! deallocate current head node 
2327 
previous => current ! make current point to new head 
2328 
END DO 
2329  
2330 
previous => Gauss_end ! make current point to head of list 
2331 
DO 
2332 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2333 
current => previous%next ! make list point to next node of head 
2334 
DEALLOCATE(previous) ! deallocate current head node 
2335 
previous => current ! make current point to new head 
2336 
END DO 
2337  
2338  
2339 
END IF 
2340  
2341 
DEALLOCATE(GeomTmp, Hess, GradTmp) 
2342 
IF (COORD.EQ.'ZMAT') DEALLOCATE(dzdc) 
2343 
IF (COORD.EQ.'BAKER') THEN 
2344 
DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t) 
2345 
DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local) 
2346 
END IF 
2347  
2348 
WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt 
2349  
2350 
Close(IOIN) 
2351 
Close(IOOUT) 
2352 
IF (AnaGeom) Close(IODAT) 
2353 
IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot) 
2354  
2355 
END PROGRAM PathOpt 