root / src / Path.f90 @ 7
Historique  Voir  Annoter  Télécharger (79,55 ko)
1 
! This programs generates and optimizes a reaction path 

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

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

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

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

629 
if (Prog.EQ.'TEST2') Prog="CHAMFRE" 
630 
if (Prog.EQ.'2D') Prog="TEST2D" 
631 
if (Prog.EQ.'TEST3') Prog="LEPS" 
632  
633  
634 
Select Case (Prog) 
635 
CASE ("EMPTY","G09","GAUSSIAN") 
636 
Prog='GAUSSIAN' 
637 
If (ProgExe=="EMPTY") ProgExe="g09" 
638 
if (CalcName=="EMPTY") CalcName="Path" 
639 
if (ISuffix=="EMPTY") ISuffix="com" 
640 
if (OSuffix=="EMPTY") OSuffix="log" 
641 
CASE ("G03") 
642 
Prog='GAUSSIAN' 
643 
If (ProgExe=="EMPTY") ProgExe="g03" 
644 
if (CalcName=="EMPTY") CalcName="Path" 
645 
if (ISuffix=="EMPTY") ISuffix="com" 
646 
if (OSuffix=="EMPTY") OSuffix="log" 
647 
CASE ("EXT") 
648 
If (ProgExe=="EMPTY") ProgExe="./Ext.exe" 
649 
if (CalcName=="EMPTY") CalcName="Ext" 
650 
if (ISuffix=="EMPTY") ISuffix="in" 
651 
if (OSuffix=="EMPTY") OSuffix="out" 
652 
CASE ('MOPAC') 
653 
If (ProgExe=="EMPTY") ProgExe="mopac" 
654 
if (CalcName=="EMPTY") CalcName="Path" 
655 
if (ISuffix=="EMPTY") ISuffix="mop" 
656 
if (OSuffix=="EMPTY") OSuffix="out" 
657 
CASE ("SIESTA") 
658 
If (ProgExe=="EMPTY") ProgExe="siesta" 
659 
if (CalcName=="EMPTY") CalcName="Path" 
660 
if (ISuffix=="EMPTY") ISuffix="fdf" 
661 
if (OSuffix=="EMPTY") OSuffix="out" 
662 
CASE ("VASP") 
663 
If (ProgExe=="EMPTY") ProgExe="VASP" 
664 
CASE ("TM","TURBOMOLE") 
665 
Prog="TURBOMOLE" 
666 
If (ProgExe=="EMPTY") ProgExe="jobex c 1 keep" 
667 
CASE ("TEST","CHAMFRE","TEST2D","LEPS") 
668 
ProgExe="" 
669 
CASE DEFAULT 
670 
WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP" 
671 
STOP 
672 
END SELECT 
673  
674 
if (Input.EQ.'EMPTY') THEN 
675 
Select Case (Prog) 
676 
CASE ("VASP") 
677 
Input=Prog 
678 
CASE ("CHAMFRE") 
679 
Input="VASP" 
680 
CASE DEFAULT 
681 
Input='XYZ' 
682 
END SELECT 
683 
WRITE(*,*) "Input has been set to the default: ",INPUT 
684 
END IF 
685  
686 
WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) 
687  
688 
IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN 
689 
WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode." 
690 
STOP 
691 
END IF 
692  
693 
if (OptGeom.GE.1) THEN 
694 
FPrintGeom=.FALSE. 
695 
If (GeomFile/='EMPTY') THEN 
696 
FPrintGeom=.TRUE. 
697 
END IF 
698 
IF (IGeomRef.LE.0) THEN 
699 
IGeomRef=OptGeom 
700 
ELSE 
701 
IF (IGeomRef.NE.OptGeom) THEN 
702 
WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
703 
WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
704  
705 
WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
706 
WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
707 
STOP 
708 
END IF 
709 
END IF 
710 
END IF 
711  
712 
ALLOCATE(Cart(Nat)) 
713 
Cart=0 
714  
715 
ALLOCATE(Frozen(Nat)) 
716 
Frozen=0 
717  
718 
IF (FCart) THEN 
719 
! We rewind the file to be sure to browse all of it to read the Cart namelist 
720 
REWIND(IOIN) 
721 
List=0 
722 
READ(IOIN,cartlist) 
723 
IF (COORD.NE.'CART') THEN 
724 
Cart=List(1:Nat) 
725 
if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart 
726 
ELSE 
727 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
728 
WRITE(*,*) "I will discard the cartlist" 
729 
Cart=0 
730 
END IF 
731 
END IF 
732  
733 
if (FFrozen) THEN 
734  
735 
REWIND(IOIN) 
736 
List=0 
737 
READ(IOIN,Frozenlist) 
738 
Frozen=List(1:Nat) 
739 
if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen 
740 
END IF 
741  
742 
If (FCart) Call Expand(Cart,NCart,Nat) 
743 
IF (FFrozen) Call Expand(Frozen,NFroz,Nat) 
744  
745 
if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart 
746 
if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen 
747  
748 
if (AnaGeom) THEN 
749 
REWIND(IOIN) 
750 
READ(IOIN,AnaList,IOSTAT=IoS) 
751 
IF (IOS==0) THEN ! we have a AnaList namelist 
752 
Call ReadAnaList 
753 
END IF 
754 
IF (AnaFile=="EMPTY") THEN 
755 
AnaFile=TRIM(PathName) // '.dat' 
756 
END IF 
757 
OPEN(IODat,File=AnaFile) 
758 
END IF 
759  
760 

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

1791 
DO Iat=max(NCart+1,4),Nat 
1792 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
1793 
OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & 
1794 
OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & 
1795 
OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi 
1796 
Idx=Idx+3 
1797 
END DO 
1798 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
1799  
1800 
CASE ('HYBRID','CART') 
1801 
XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
1802 
CASE DEFAULT 
1803 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP' 
1804 
STOP 
1805 
END SELECT 
1806  
1807 
IF (debug) THEN 
1808 
WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt 
1809 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1810 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1811 
ELSE 
1812 
DO Iat=1,Nat 
1813 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1814 
XyzGeomF(IGeom,1:3,Iat) 
1815 
END DO 
1816 
END IF 
1817 
END IF ! matches if (debug) THEN 
1818  
1819 
! We project out the tangential part of the gradient to know if we are done 
1820 
GradTmp=Grad(IGeom,:) 
1821 
NormGrad=sqrt(dot_product(GradTmp,GradTmp)) 
1822 
if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad 
1823 
SELECT CASE (COORD) 
1824 
CASE ('ZMAT','MIXED','BAKER') 
1825 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
1826 
Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:)) 
1827 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
1828 
Ca=S/(sqrt(Norm)*NormGrad) 
1829 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
1830 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
1831 
CASE ('CART','HYBRID') 
1832 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
1833 
Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:)) 
1834 
Ca=S/(sqrt(Norm)*NormGrad) 
1835 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
1836 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
1837 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
1838 
CASE DEFAULT 
1839 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP' 
1840 
STOP 
1841 
END SELECT 
1842 
IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi 
1843 
NormGrad=sqrt(DOT_Product(GradTmp,GradTmp)) 
1844 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
1845 
IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini 
1846  
1847 
END DO ! matches DO IGeom=2,NGeomF1 
1848 
! we save the old gradients 
1849 
GradOld=Grad 
1850 
EPathOld=Energies 
1851  
1852  
1853  
1854 
! Shall we reparameterize the path ? 
1855 
! PFL 2007/Apr/10 > 
1856 
! We call PathCreate in all cases, it will take care of the 
1857 
! reparameterization, as well as calculating the tangents. 
1858 
! IF (MOD(IOpt,IReparam).EQ.0) THEN 
1859 
XyzGeomI=XyzGeomF 
1860 
IntCoordI=IntCoordF 
1861  
1862 
Call PathCreate(IOpt) 
1863 
! END IF 
1864 
! < PFL 2007/Apr/10 
1865  
1866 
if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)') 'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini 
1867 
IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN 
1868  
1869 
! We have the new path, we calculate its energy and gradient 
1870  
1871 
Call EgradPath(IOpt) 
1872 
!IF(IOPT .EQ. 10) Then 
1873 
! Print *, 'Stopping at my checkpoint.' 
1874 
! STOP !This is my temporary checkpoint. 
1875 
!ENDIF 
1876  
1877 
! PFL 31 Mar 2008 > 
1878 
! Taken from v3.99 but we added a flag... 
1879 
! Added in v3.99 : MaxStep is modified according to the evolution of energies 
1880 
! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
1881 
! else it is decreased by 0.8 
1882  
1883 
if (DynMaxStep) THEN 
1884 
if (debug) WRITE(*,*) "Dynamically updating the Maximum step" 
1885 
if (OptReac) THEN 
1886 
If (Energies(1)<EPathOld(1)) Then 
1887 
MaxStep(1)=min(1.1*MaxStep(1),2.*SMax) 
1888 
ELSE 
1889 
MaxStep(1)=max(0.8*MaxStep(1),SMax/2.) 
1890 
END IF 
1891 
END IF 
1892  
1893 
if (OptProd) THEN 
1894 
If (Energies(NGeomF)<EPathOld(NGeomF)) Then 
1895 
MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax) 
1896 
ELSE 
1897 
MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.) 
1898 
END IF 
1899 
END IF 
1900  
1901 
DO IGeom=2,NGeomF1 
1902 
If (Energies(IGeom)<EPathOld(IGeom)) Then 
1903 
MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax) 
1904 
ELSE 
1905 
MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.) 
1906 
END IF 
1907 
END DO 
1908 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
1909 
END IF ! test on DynMaxStep 
1910 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
1911 
! < PFL 31 MAr 2008 
1912 
! Also XyzGeomF is updated in EgradPath for Baker Case. 
1913 
! It should get updated for other cases too. 
1914  
1915 
DO IGeom=1,NGeomF 
1916 
SELECT CASE (COORD) 
1917 
CASE('ZMAT') 
1918 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
1919 
GeomTmp=IntCoordF(IGeom,:) 
1920 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) 
1921 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
1922 
OrderInv(indzmat(2,2)),Geomtmp(1) 
1923 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
1924 
OrderInv(indzmat(3,2)),Geomtmp(2), & 
1925 
OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi 
1926 
Idx=4 
1927 
DO Iat=4,Nat 
1928 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
1929 
OrderInv(indzmat(iat,2)),Geomtmp(Idx), & 
1930 
OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, & 
1931 
OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi 
1932 
Idx=Idx+3 
1933 
END DO 
1934 
CASE('BAKER') 
1935 
!WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
1936 
!GeomTmp=IntCoordF(IGeom,:) 
1937 
CASE('CART','HYBRID','MIXED') 
1938 
WRITE(IOOUT,'(1X,I5)') Nat 
1939 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
1940 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1941 
DO I=1,Nat 
1942 
Iat=I 
1943 
If (renum) Iat=OrderInv(I) 
1944 
WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat2:3*Iat) 
1945 
END DO 
1946 
CASE DEFAULT 
1947 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP' 
1948 
STOP 
1949 
END SELECT 
1950 
END DO ! matches DO IGeom=1,NGeomF 
1951  
1952 
Call Write_path(IOpt) 
1953 
! Shall we analyze the geometries ? 
1954 
IF (AnaGeom) Call AnaPath(Iopt,IoDat) 
1955  
1956 
! We update the Hessian matrices 
1957 
IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord 
1958 
! First using the displacement from the old path 
1959 
IGeom0=2 
1960 
IGeomF=NGeomF1 
1961 
IF (OptReac) IGeom0=1 
1962 
IF (OptProd) IGeomF=NGeomF 
1963  
1964 
IF (FCalcHess) THEN 
1965 
DO IGeom=IGeom0,IGeomF 
1966 
SELECT CASE (COORD) 
1967 
CASE ('ZMAT','MIXED','BAKER') 
1968 
GeomTmp2=IntCoordF(IGeom,:) 
1969 
CASE ('CART','HYBRID') 
1970 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1971 
CASE DEFAULT 
1972 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
1973 
STOP 
1974 
END SELECT 
1975 
Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt) 
1976 
END DO 
1977 
ELSE 
1978 
DO IGeom=IGeom0,IGeomF 
1979 
GeomTmp=GeomOld_all(IGeom,:) 
1980 
SELECT CASE (COORD) 
1981 
CASE ('ZMAT','MIXED','BAKER') 
1982 
GeomTmp2=IntCoordF(IGeom,:) 
1983 
CASE ('CART','HYBRID') 
1984 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1985 
CASE DEFAULT 
1986 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
1987 
STOP 
1988 
END SELECT 
1989 
GeomTmp=GeomTmp2GeomTmp 
1990 
GradTmp=Grad(IGeom,:)GradOld(IGeom,:) 
1991 
HessTmp=Hess(IGeom,:,:) 
1992 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
1993 
Hess(IGeom,:,:)=HessTmp 
1994 
END DO ! matches DO IGeom=IGeom0,IGeomF 
1995  
1996 
! Second using the neighbour points 
1997 
IF (HupNeighbour) THEN 
1998 
! Only one point for end points. 
1999 
IF (OptReac) THEN 
2000 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2001 
GeomTmp=IntCoordF(1,:)IntCoordF(2,:) 
2002 
ELSE 
2003 
GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))Reshape(XyzGeomF(2,:,:),(/3*Nat/)) 
2004 
END IF 
2005 
GradTmp=Grad(1,:)Grad(2,:) 
2006 
HessTmp=Hess(1,:,:) 
2007 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2008 
Hess(1,:,:)=HessTmp 
2009 
END IF 
2010  
2011 
IF (OptProd) THEN 
2012 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2013 
GeomTmp=IntCoordF(NGeomF,:)IntCoordF(NGeomF1,:) 
2014 
ELSE 
2015 
GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))Reshape(XyzGeomF(NGeomF1,:,:),(/3*Nat/)) 
2016 
END IF 
2017 
GradTmp=Grad(NGeomF,:)Grad(NGeomF1,:) 
2018 
HessTmp=Hess(NGeomF,:,:) 
2019 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2020 
Hess(NGeomF,:,:)=HessTmp 
2021 
END IF 
2022  
2023 
! Two points for the rest of the path 
2024 
DO IGeom=2,NGeomF1 
2025 
IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN 
2026 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom+1,:) 
2027 
ELSE 
2028 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) 
2029 
END IF 
2030 
GradTmp=Grad(IGeom,:)Grad(IGeom+1,:) 
2031 
HessTmp=Hess(IGeom,:,:) 
2032 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2033  
2034 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2035 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom1,:) 
2036 
ELSE 
2037 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom1,:,:),(/3*Nat/)) 
2038 
END IF 
2039 
GradTmp=Grad(IGeom,:)Grad(IGeom1,:) 
2040  
2041 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2042 
Hess(IGeom,:,:)=HessTmp 
2043 
END DO ! matches DO IGeom=2,NGeomF1 
2044 
END IF !matches IF (HupNeighbour) THEN 
2045 
END IF ! matches IF (FCalcHess) 
2046 
END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc) 
2047 
END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) 
2048  
2049 
! IF (PROG=="GAUSSIAN") THEN 
2050 
! DEALLOCATE(Gauss_paste) 
2051 
! DEALLOCATE(Gauss_root) 
2052 
! DEALLOCATE(Gauss_comment) 
2053 
! DEALLOCATE(Gauss_end) 
2054 
! END IF 
2055  
2056 
DEALLOCATE(XyzGeomF, IntCoordF) 
2057 
DEALLOCATE(XyzGeomI, IntCoordI) 
2058 
DEALLOCATE(XyzTangent,IntTangent) 
2059 
DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt) 
2060 
DEALLOCATE(GeomOld_all) 
2061  
2062 
if (PROG=="GAUSSIAN") THEN 
2063 
! We deallocate the Gauss_XX lists 
2064 
! transverse the list and deallocate each node 
2065 
previous => Gauss_root ! make current point to head of list 
2066 
DO 
2067 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2068 
current => previous%next ! make list point to next node of head 
2069 
DEALLOCATE(previous) ! deallocate current head node 
2070 
previous => current ! make current point to new head 
2071 
END DO 
2072  
2073 
previous => Gauss_comment ! make current point to head of list 
2074 
DO 
2075 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2076 
current => previous%next ! make list point to next node of head 
2077 
DEALLOCATE(previous) ! deallocate current head node 
2078 
previous => current ! make current point to new head 
2079 
END DO 
2080  
2081 
previous => Gauss_end ! make current point to head of list 
2082 
DO 
2083 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2084 
current => previous%next ! make list point to next node of head 
2085 
DEALLOCATE(previous) ! deallocate current head node 
2086 
previous => current ! make current point to new head 
2087 
END DO 
2088  
2089  
2090 
END IF 
2091  
2092 
DEALLOCATE(GeomTmp, Hess, GradTmp) 
2093 
IF (COORD.EQ.'ZMAT') DEALLOCATE(dzdc) 
2094 
IF (COORD.EQ.'BAKER') THEN 
2095 
DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t) 
2096 
DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local) 
2097 
END IF 
2098  
2099 
WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt 
2100  
2101 
Close(IOIN) 
2102 
Close(IOOUT) 
2103 
IF (AnaGeom) Close(IODAT) 
2104  
2105 
END PROGRAM PathOpt 