root / src / Path.f90 @ 10
History  View  Annotate  Download (89.9 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,TChk 
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 
FPBC=.TRUE. 
835 
IPer=3 
836 
Kabeg=1 
837 
kaEnd=1 
838 
kbBeg=1 
839 
kbEnd=1 
840 
kcBeg=1 
841 
kcEnd=1 
842 
CASE ("TM","TURBOMOLE") 
843 
Prog="TURBOMOLE" 
844 
If (ProgExe=="EMPTY") ProgExe="jobex c 1 keep" 
845 
UnitProg="au" 
846 
ConvE=au2kcal 
847 
FPBC=.FALSE. 
848 
CASE ("TEST","CHAMFRE","TEST2D","LEPS") 
849 
ProgExe="" 
850 
UnitProg="au" 
851 
ConvE=au2kcal 
852 
CASE DEFAULT 
853 
WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP" 
854 
STOP 
855 
END SELECT 
856  
857 
if (Input.EQ.'EMPTY') THEN 
858 
Select Case (Prog) 
859 
CASE ("VASP") 
860 
Input=Prog 
861 
CASE ("CHAMFRE") 
862 
Input="VASP" 
863 
CASE DEFAULT 
864 
Input='XYZ' 
865 
END SELECT 
866 
WRITE(*,*) "Input has been set to the default: ",INPUT 
867 
END IF 
868  
869 
WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) 
870  
871 
IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN 
872 
WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode." 
873 
STOP 
874 
END IF 
875  
876 
if (OptGeom.GE.1) THEN 
877 
FPrintGeom=.FALSE. 
878 
If (GeomFile/='EMPTY') THEN 
879 
FPrintGeom=.TRUE. 
880 
END IF 
881 
IF (IGeomRef.LE.0) THEN 
882 
IGeomRef=OptGeom 
883 
ELSE 
884 
IF (IGeomRef.NE.OptGeom) THEN 
885 
WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
886 
WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
887  
888 
WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
889 
WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
890 
STOP 
891 
END IF 
892 
END IF 
893 
END IF 
894  
895 
ALLOCATE(Cart(Nat)) 
896 
Cart=0 
897  
898 
ALLOCATE(Frozen(Nat)) 
899 
Frozen=0 
900  
901 
IF (FCart) THEN 
902 
! We rewind the file to be sure to browse all of it to read the Cart namelist 
903 
REWIND(IOIN) 
904 
List=0 
905 
READ(IOIN,cartlist) 
906 
IF (COORD.NE.'CART') THEN 
907 
Cart=List(1:Nat) 
908 
if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart 
909 
ELSE 
910 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
911 
WRITE(*,*) "I will discard the cartlist" 
912 
Cart=0 
913 
END IF 
914 
END IF 
915  
916 
if (FFrozen) THEN 
917  
918 
REWIND(IOIN) 
919 
List=0 
920 
READ(IOIN,Frozenlist) 
921 
Frozen=List(1:Nat) 
922 
if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen 
923 
END IF 
924  
925 
If (FCart) Call Expand(Cart,NCart,Nat) 
926 
IF (FFrozen) Call Expand(Frozen,NFroz,Nat) 
927  
928 
if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart 
929 
if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen 
930  
931 
if (AnaGeom) THEN 
932 
REWIND(IOIN) 
933 
READ(IOIN,AnaList,IOSTAT=IoS) 
934 
IF (IOS==0) THEN ! we have a AnaList namelist 
935 
Call ReadAnaList 
936 
ELSE 
937 
! No AnaList namelist, we will print only the energy 
938 
FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)' 
939 
END IF 
940 
IF (AnaFile=="EMPTY") THEN 
941 
AnaFile=TRIM(PathName) // '.dat' 
942 
END IF 
943 
OPEN(IODat,File=AnaFile) 
944 
IF (GplotFile=="EMPTY") THEN 
945 
GplotFile=TRIM(PathName) // '.gplot' 
946 
END IF 
947 

948 
OPEN(IODat,File=AnaFile) 
949 
If (OptGeom<=0) THEN 
950 
OPEN(IOGplot,File=GPlotFile) 
951 
WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot persist" 
952 
WRITE(IoGplot,'(A)') " set pointsize 2" 
953 
WRITE(IoGplot,'(A)') " xcol=1" 
954 
WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3 
955 
END IF 
956 
END IF 
957  
958 

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

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