root / src / Path.f90 @ 8
History  View  Annotate  Download (89.2 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 namelis. 
107 
! v3.96 
108 
! I have added a test for the geometry step: it does not allow valence angle to 
109 
! be negative or bigger than 180. I hope that this might help geometry optimization 
110 
! when dealing with nearly linear molecules. 
111 
! v3.961 
112 
! The options cart and frozen are now LOGICAL: 
113 
! Fcart=T indicates that a &cartlist namelist should be read ! 
114 
! Ffrozen=T indicates that a &frozenlist namelist should be read ! 
115 
! TO DO: Autocart=T indicates that the atoms described in cartesian coordiantes should 
116 
! be found automatically 
117 
! v3.97 
118 
! Baker has been implemented by P. Dayal: we are interpolating Xint for now. 
119 
! TO do: interpolate Umat ? 
120 
!  
121 
! v3.971 
122 
! AutoCart=T is being implemented. Goal: having a 'black box' version for Vasp users. 
123 
! Vmd: True if user wants to watch the path using VMD. Only used for VASP for now. 
124 
! 
125 
! v3.98 
126 
! Calc_zmat has been replaced by Calc_zmat_frag which is more robust. 
127 
! Still missing: linear molecules and some tests. 
128 
! 
129 
! v3.981 
130 
! * Added a non document flag in .valid file that switch on the numerical 
131 
! calculation of the hessian instead of the update. 
132 
! * Added HesUpd variable in the namelist. Default value is MS for path 
133 
! optimization because it does not imposes a positive hessian, and BFGS 
134 
! for geometry optimization 
135 
! v 4.0 
136 
! * This version has been made public within the CARTE project. 
137 
! * Added: 
138 
!  Step_method to choose the optimization method: RFO or Diis 
139 
!  DynMaxStep: if TRUE, the maximum size of a step is updated at each step. 
140 
! If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. 
141 
! It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0] 
142 
! 
143 
!!!!!!!!!!!!!!!!!!!!!! 
144 
! v4.1 
145 
! * Para mode has been partly implemented for Gaussian. 
146 
! VaspMode is thus now 'RunMode' and can be Serial or Para 
147 
! Only Gaussian and VASP are supported for Para mode now. 
148 
! v4.12 
149 
! Some bugs in unconstrained zmat construction are corrected. 
150 
! v4.13 
151 
! Prakash has implemented the GEDIIS optimization procedure. 
152 
! v4.14 
153 
! There is a small problem for some inputs in VASP, where the 
154 
! last geometry is not the one given by the user. 
155 
! It seems to come from the fact that some people try to freeze 
156 
! some unfrozen atoms 
157 
! So some tests have been added to check that frozen atoms do not move. 
158 
!!! 
159 
! v4.15 
160 
! There were some bugs when reading and expanding cartlist and frozenlist 
161 
! I have changed the way frozen atoms are compared by adding partial alignment. 
162 
! v4.16 
163 
! Some bugs were corrected. 
164 
! v4.17 
165 
! Added MOPAC as a program. 
166 
! v4.175 
167 
! Added Three more parameters for defining the input and output names for 
168 
! the calculations. 
169 
! CalcName: Prefix for the files used for the energy and gradient calculations 
170 
! ISuffix: Suffix for the input file 
171 
! OSuffix: suffix for the output file. 
172 
! This maters for Gaussian for example. 
173 
! 
174 
! v4.177 (P) 2010 Nov 22 
175 
!  We add TurboMole as a new external code. 
176 
!  Subtle change for Input: the default is XYZ whatever the code, except 
177 
! for VASP. 
178 
! 
179 
! v4.178 (P) 2010 Nov 28 
180 
!  We add the possibility to call CHAMFRE reactive force field using 
181 
! prog=TEST2 or CHAMFRE 
182 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
183 
! OpenPath v1.4 2011 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 
PROGRAM PathOpt 
209  
210 
use VarTypes 
211 
use Path_module 
212 
use Io_module 
213  
214 
IMPLICIT NONE 
215  
216 
CHARACTER(132) :: FileIn, FileOut, Line,Line2 
217 
CHARACTER(32) :: Input,poscar 
218  
219 
INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF 
220 
INTEGER(KINT) :: IOpt 
221 
INTEGER(KINT) :: Idx, LineL,LenLine 
222 
INTEGER(KINT) :: N3at, NFree, Nfr 
223 
! Temporary integer counter 
224 
INTEGER(KINT) :: NTmp 
225  
226 
INTEGER(KINT) :: List(2*MaxFroz) 
227  
228 
REAL(KREAL) :: E,FactStep,B(3) !E=Energy 
229 
REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad 
230 
REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E 
231 
REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord 
232 
REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord) 
233 
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord 
234 
REAL(KREAL), ALLOCATABLE :: HessTmp(:,:) 
235 
REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim) 
236 
REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat6)) 
237 
LOGICAL, SAVE :: allocation_flag 
238  
239 
! Flag to see if frozen atoms are frozen (see below) 
240 
LOGICAL :: FChkFrozen 
241  
242 
! Energies for all points along the path at the previous iteration 
243 
REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF, where it is deallocated: Prakash 
244 
! Maximum step allowed for this geometry 
245 
REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF, where it is deallocated: Prakash 
246  
247 
! these are used to read temporary coordinates 
248 
REAL(KREAL) :: xtmp,ytmp,ztmp 
249  
250 
LOGICAL :: FFrozen,FCart 
251  
252 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
253 
! 
254 
! To analyse the frozen atoms 
255 
! 
256 
REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat 
257 
REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat 
258 
REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist 
259 
LOGICAL, ALLOCATABLE :: ListAtFroz(:) 
260 
INTEGER(KINT) :: Iat1, Iat2, Iat3 
261  
262 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
263 
! 
264 
! For VASP 
265 
! 
266 
INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na 
267 
INTEGER(KINT) :: NbType, NbTypeUser 
268  
269 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
270  
271 
LOGICAL :: TChk 
272 
LOGICAL :: Debug, Fini,Flag_Opt_Geom 
273  
274 
! INTEGER(KINT), EXTERNAL :: Iargc 
275 
INTEGER(KINT), EXTERNAL :: ConvertNumAt 
276  
277 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
278 
! 
279 
! For Debugging purposes 
280 
! 
281 
LOGICAL :: FCalcHess 
282 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
283  
284  
285 
INTERFACE 
286 
function valid(string) result (isValid) 
287 
CHARACTER(*), intent(in) :: string 
288 
logical :: isValid 
289 
END function VALID 
290  
291 
SUBROUTINE Read_Geom(Input) 
292 
CHARACTER(32), INTENT(IN) :: input 
293 
END SUBROUTINE READ_GEOM 
294  
295 
subroutine hupdate_all (n, dx, dgrad, HessOld) 
296  
297 
IMPLICIT NONE 
298  
299 
INTEGER, PARAMETER :: KINT=KIND(1) 
300 
INTEGER, PARAMETER :: KREAL=KIND(1.0D0) 
301  
302 
INTEGER(KINT), INTENT(IN) :: n 
303 
real(KREAL) :: dx(n), dgrad(n), HessOld(n*n) 
304 
END subroutine hupdate_all 
305  
306 
SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess) 
307  
308 
IMPLICIT NONE 
309 

310 
INTEGER, PARAMETER :: KINT=KIND(1) 
311 
INTEGER, PARAMETER :: KREAL=KIND(1.D0) 
312 

313 
INTEGER(KINT), INTENT(IN) :: NFREE 
314 
REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree) 
315 
REAL(KREAL), INTENT(IN) :: ds(Nfree) 
316 
REAL(KREAL), INTENT(IN) :: dGrad(Nfree) 
317 

318 
END subroutine Hinvup_BFGS_new 
319  
320  
321 
END INTERFACE 
322  
323 
Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, & 
324 
MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, & 
325 
SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, & 
326 
input,prog, ProgExe,RunMode, AtTypes,poscar,PathName, & 
327 
OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, & 
328 
FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, & 
329 
NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, & 
330 
NGintMax, Align, CalcEReac,CalcEProd 
331  
332 
Namelist/cartlist/list 
333 
Namelist/frozenlist/list 
334  
335  
336 
Flag_Opt_Geom = .FALSE. ! Initialized. 
337  
338 
!!!!!!!!!!!!!!!!!!!!!!!!! 
339 
! 
340 
! We print the Version info in file 
341 
! 
342 
WRITE(*,'(1X,A)') "OpenPath v1.4 (c) PFL/PD 20072011" 
343  
344 
! We read the files names 
345 
! SELECT CASE (iargc()) 
346 
SELECT CASE (command_argument_count()) 
347 
CASE (2) 
348 
call getarg(1,FileIn) 
349 
OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') 
350 
call getarg(2,FileOut) 
351 
OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN') 
352 
CASE (1) 
353 
call getarg(1,FileIn) 
354 
IF ((FileIn=="help").OR.(FileIn=="help")) THEN 
355 
WRITE(*,*) "Path minihelp" 
356 
WRITE(*,*) "" 
357 
WRITE(*,*) "" 
358 
WRITE(*,*) "Use: Path Input_file Output_file" 
359 
WRITE(*,*) "Input_file starts with a Namelist called path" 
360 
WRITE(*,*) "" 
361 
WRITE(*,*) "Compulsory variables are:" 
362 
WRITE(*,*) "NGeomi: Number of geometries defining the Initial path" 
363 
WRITE(*,*) "NGeomf: Number of geometries defining the Final path" 
364 
WRITE(*,*) "Nat : Number of atoms" 
365 
WRITE(*,*) "" 
366 
WRITE(*,*) "" 
367 
WRITE(*,*) "Other options" 
368 
WRITE(*,*) "Input: string that indicates the type of the input geometries." 
369 
WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp" 
370 
WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations." 
371 
WRITE(*,*) " Accepted values are: Gaussian, Vasp, Mopac or Ext" 
372 
WRITE(*,*) " In case of a Gaussian or Mopac calculations, input must be set to Cart." 
373 
WRITE(*,*) " One example of a gaussian/mopac input should be added at the end of the" 
374 
WRITE(*,*) " input file.See example file Test_G03.path or Test_Mopac.path" 
375 
WRITE(*,*) " In the case of a VASP calculation, if input is set to Cart, then" 
376 
WRITE(*,*) " the preamble of a VASP calculation should be added at the end of" 
377 
WRITE(*,*) " the input file. See example file Test_VASP_cart.path" 
378 
WRITE(*,*) " In the case of a VASP calculation, one should also give value of the " 
379 
WRITE(*,*) " Runmode variable" 
380 
WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy" 
381 
WRITE(*,*) " and gradient of the whole path or not. If one wants to use VASP," 
382 
WRITE(*,*) " Runmode must be set to PARA." 
383 
WRITE(*,*) "PathName: Prefix used to save the path. Default is Path" 
384 
WRITE(*,*) "Poscar: string that will be used as the prefix for the different " 
385 
WRITE(*,*) " POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.," 
386 
WRITE(*,*) " not used for internal calculations." 
387 
WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations" 
388 
WRITE(*,*) " ISuffix: Suffix for the input file." 
389 
WRITE(*,*) " OSuffix: suffix for the output file." 
390 
WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates." 
391 
WRITE(*,*) " Valid only for Coord=Zmat, Hybrid or Mixed" 
392 
WRITE(*,*) "Fact: REAL used to define if two atoms are linked." 
393 
WRITE(*,*) " If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked." 
394 
WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info." 
395 
WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:" 
396 
WRITE(*,*) "  CART (or Xyz): works in cartesian" 
397 
WRITE(*,*) "  Zmat: works in internal coordinates (Zmat)" 
398 
WRITE(*,*) "  Mixed: frozen atoms, as well as atoms defined by the " 
399 
WRITE(*,*) " 'cart' array(see below) are describe in CARTESIAN, whereas the" 
400 
WRITE(*,*) " others are described in Zmat" 
401 
WRITE(*,*) "  Baker: use of Baker coordinates, also called delocalized internal coordinates" 
402 
WRITE(*,*) "  Hybrid: geometries are described in zmat, but the gradient are used in cartesian" 
403 
WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
404 
WRITE(*,*) "  RFO: Rational function optimization" 
405 
WRITE(*,*) "  GDIIS: Geometry optimization using direct inversion in the iterative supspace" 
406 
! WRITE(*,*) " HesUpd: Deprecated. Use Hupdate." 
407 
WRITE(*,*) " HUpdate: method to update the hessian. By default, it is MurtaghSargent" 
408 
WRITE(*,*) " Except for geometry optimization where it is BFGS." 
409 
WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization" 
410 
WRITE(*,*) "Smax: Maximum length of a step during path optimization" 
411 
WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary" 
412 
WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken" 
413 
WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path" 
414 
WRITE(*,*) " that is kept. FTan=1. corresponds to the full displacement, " 
415 
WRITE(*,*) " whereas FTan=0. gives a displacement orthogonal to the path." 
416 
WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization." 
417 
WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path." 
418 
WRITE(*,*) " This option indicates the first step where spline interpolation is used." 
419 
WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, " 
420 
WRITE(*,*) " it might happen that an atom moves out of the unit cell." 
421 
WRITE(*,*) " Path detects this by comparing the displacement to boxtol:" 
422 
WRITE(*,*) " if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5" 
423 
WRITE(*,*) "" 
424 
WRITE(*,*) "Arrays:" 
425 
WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements." 
426 
WRITE(*,*) " You can modify it using, rcov(6)=0.8." 
427 
WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements." 
428 
WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations." 
429 
WRITe(*,*) "If not given, Path will read the POTCAR file." 
430 
WRITE(*,*) "" 
431 
WRITE(*,*) "Flags:" 
432 
WRITE(*,*) "MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE." 
433 
WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time." 
434 
WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products." 
435 
WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants." 
436 
WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para" 
437 
WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para" 
438 
WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
439 
WRITE(*,*) "Hinv: if True, then Hessian inversed is used." 
440 
WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations." 
441 
WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path." 
442 
WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
443 
WRITE(*,*) " If True, a &frozenlist namelist containing the list of frozen atoms must be given." 
444 
WRITE(*,*) " If VASP is used, and frozen is not given" 
445 
WRITE(*,*) " here, the program will use the F flags of the input geometry" 
446 
WRITE(*,*) "FCart: True if one wants to describe some atoms using cartesian coordinates. " 
447 
WRITE(*,*) " *** Only used in 'mixed' calculations. ***" 
448 
WRITE(*,*) " If True, a &cartlist namelist containing the list of cart atoms must be given." 
449 
WRITE(*,*) " By default, only frozen atoms are described in cartesian coordinates." 
450 
WRITE(*,*) "" 
451 
WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry." 
452 
WRITE(*,*) " If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. " 
453 
WRITE(*,*) " It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]" 
454 
WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms." 
455 
WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now" 
456 
WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files." 
457 
WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP." 
458 
WRITE(*,*) "" 
459 
STOP 
460 
END IF 
461 
OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') 
462 
IOOUT=6 
463 
CASE (0) 
464 
IOIN=5 
465 
IOOUT=6 
466 
END SELECT 
467  
468  
469  
470 
! We initiliaze variables 
471 
Pi=dacos(1.0d0) 
472 
Nat=1 
473 
Fact=1.1 
474 
IGeomRef=1 
475 
NGeomI=1 
476 
NGeomF=8 
477 
IReparam=5 
478 
IReparamT=1 
479 
ISpline=5 
480 
debug=.False. 
481 
MW=.TRUE. 
482 
bohr=.False. 
483 
Coord='ZMAT' 
484 
Step_method='RFO' 
485 
DebugFile='Path.valid' 
486 
PathName="Path" 
487 
MaxCyc=50 
488 
SMax=0.1D0 
489 
SThresh=1. 
490 
GThresh=1. 
491 
FTan=0. 
492 
Hinv=.TRUE. 
493 
IniHUp=.FALSE. 
494 
HupNeighbour=.FALSE. 
495 
HesUpd="EMPTY" 
496 
HUpdate="MS" 
497 
FFrozen=.FALSE. 
498 
FCart=.FALSE. 
499 
List=0 
500 
renum=.TRUE. 
501 
OptReac=.FALSE. 
502 
OptProd=.FALSE. 
503 
CalcEReac=.FALSE. 
504 
CalcEProd=.FALSE. 
505 
EReac=0.d0 
506 
EProd=0.d0 
507 
OptGeom=1 
508 
PathOnly=.FALSE. 
509 
Prog='EMPTY' 
510 
ProgExe='EMPTY' 
511 
Runmode='SERIAL' 
512 
Input='EMPTY' 
513 
Poscar="POSCAR" 
514 
AutoCart=.FALSE. 
515 
VMD=.TRUE. 
516 
WriteVASP=.FALSE. 
517 
DynMaxStep=.FALSE. 
518 
CalcName="EMPTY" 
519 
ISuffix="EMPTY" 
520 
OSuffix="EMPTY" 
521 
NGintMax=10 
522 
Align=.TRUE. 
523  
524 
! Number of point used to compute the length of the path, 
525 
! and to reparameterize the path 
526 
NMaxPtPath=1 
527 
FrozTol=1e4 
528  
529 
! Read the namelist 
530 
read (IOIN,path) 
531  
532 
debug=valid("Main") 
533 
FCalcHess=valid("CalcHess") 
534 
PathName=AdjustL(PathName) 
535 
Coord=AdjustL(Coord) 
536 
CALL UpCase(coord) 
537  
538 
Step_method=AdjustL(Step_method) 
539 
CALL UpCase(Step_method) 
540  
541 
Prog=AdjustL(Prog) 
542 
CALL UpCase(prog) 
543  
544 
Input=AdjustL(Input) 
545 
CALL UpCase(input) 
546  
547 
Runmode=AdjustL(Runmode) 
548 
Call UpCase(Runmode) 
549 
IF (Runmode(1:4).EQ.'PARA') Runmode="PARA" 
550  
551 
if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN 
552 
WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'" 
553 
WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE" 
554 
CalcEReac=.FALSE. 
555 
CalcEProd=.FALSE. 
556 
END IF 
557  
558 
If (IReparamT<=0) IReparamT=IReparam 
559  
560 
! We take care of HesUpd flag 
561 
! this is needed because some users might still use it 
562 
if (HesUpd.NE."EMPTY") THEN 
563 
! We are pityless: 
564 
WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead" 
565 
WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead" 
566 
STOP 
567 
END IF 
568  
569  
570 
IF (COORD.EQ.'CART') THEN 
571 
Renum=.FALSE. 
572 
IGeomRef=1 
573 
END IF 
574  
575 
if (Prog.EQ.'EMPTY') THEN 
576 
Prog='GAUSSIAN' 
577 
If (ProgExe=="EMPTY") ProgExe="g09" 
578 
END IF 
579 
if (Prog.EQ.'G03') THEN 
580 
Prog='GAUSSIAN' 
581 
If (ProgExe=="EMPTY") ProgExe="g03" 
582 
END IF 
583 
if (Prog.EQ.'G09') THEN 
584 
Prog='GAUSSIAN' 
585 
If (ProgExe=="EMPTY") ProgExe="g09" 
586 
END IF 
587  
588 
if (Prog.EQ.'TM') Prog="TURBOMOLE" 
589 
if (Prog.EQ.'TEST2') Prog="CHAMFRE" 
590  
591  
592 
Select Case (Prog) 
593 
CASE ("EXT") 
594 
if (CalcName=="EMPTY") CalcName="Ext" 
595 
if (ISuffix=="EMPTY") ISuffix="in" 
596 
if (OSuffix=="EMPTY") OSuffix="out" 
597 
CASE ('MOPAC') 
598 
if (CalcName=="EMPTY") CalcName="Path" 
599 
if (ISuffix=="EMPTY") ISuffix="mop" 
600 
if (OSuffix=="EMPTY") OSuffix="out" 
601 
CASE ("GAUSSIAN") 
602 
if (CalcName=="EMPTY") CalcName="Path" 
603 
if (ISuffix=="EMPTY") ISuffix="com" 
604 
if (OSuffix=="EMPTY") OSuffix="log" 
605 
CASE DEFAULT 
606 
if (CalcName=="EMPTY") CalcName="Path" 
607 
if (ISuffix=="EMPTY") ISuffix="in" 
608 
if (OSuffix=="EMPTY") OSuffix="out" 
609 
END SELECT 
610  
611 
if (Input.EQ.'EMPTY') THEN 
612 
Select Case (Prog) 
613 
CASE ("VASP") 
614 
Input=Prog 
615 
CASE ("CHAMFRE") 
616 
Input="VASP" 
617 
CASE DEFAULT 
618 
Input='XYZ' 
619 
END SELECT 
620 
WRITE(*,*) "Input has been set to the default: ",INPUT 
621 
END IF 
622  
623 
WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) 
624  
625 
IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN 
626 
WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode." 
627 
STOP 
628 
END IF 
629  
630 
if (OptGeom.GE.1) THEN 
631 
HUpdate="BFGS" 
632 
IF (IGeomRef.LE.0) THEN 
633 
IGeomRef=OptGeom 
634 
ELSE 
635 
IF (IGeomRef.NE.OptGeom) THEN 
636 
WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
637 
WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
638  
639 
WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !" 
640 
WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom 
641 
STOP 
642 
END IF 
643 
END IF 
644 
END IF 
645  
646 
IF (ProgExe.EQ.'EMPTY') THEN 
647 
SELECT CASE (PROG) 
648 
CASE ("GAUSSIAN") 
649 
ProgExe="g03" 
650 
CASE ("MOPAC") 
651 
ProgExe="mopac" 
652 
CASE ("EXT") 
653 
ProgExe="./Ext.exe" 
654 
CASE ("VASP") 
655 
ProgExe='VASP' 
656 
CASE ("TURBOMOLE") 
657 
ProgExe='jobex c 1 keep' 
658 
CASE ("TEST") 
659 
ProgExe="" 
660 
CASE ("CHAMFRE") 
661 
ProgExe="" 
662 
CASE DEFAULT 
663 
WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP" 
664 
STOP 
665 
END SELECT 
666 
END IF 
667  
668 
if (FCart.AND.FFrozen) THEN 
669 
! We have to be carreful because user might use any order 
670 
! to give the two lists CartList and FrozenList 
671 
READ(IOIN,'(A)') Line 
672 
Call Upcase(Line) 
673 
if (Index(Line,"CART").NE.0) THEN 
674 
ALLOCATE(Cart(Nat)) 
675 
BackSpace(IOIN) 
676 
List=0 
677 
READ(IOIN,cartlist) 
678 
IF (COORD.NE.'CART') THEN 
679 
Cart=List(1:Nat) 
680 
if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart 
681 
ELSE 
682 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
683 
WRITE(*,*) "I will discard the cartlist" 
684 
Cart=0 
685 
END IF 
686  
687 
ALLOCATE(Frozen(Nat)) 
688 
BACKSPACE(IOIN) 
689 
List=0 
690 
READ(IOIN,Frozenlist) 
691 
Frozen=List(1:Nat) 
692 
if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen 
693 
ELSE 
694 
ALLOCATE(Frozen(Nat)) 
695 
BACKSPACE(IOIN) 
696 
List=0 
697 
READ(IOIN,Frozenlist) 
698 
Frozen=List(1:Nat) 
699 
if (debug) WRITE(*,*) "DBG Path L523 Frozen",Frozen 
700 
ALLOCATE(Cart(Nat)) 
701 
BackSpace(IOIN) 
702 
List=0 
703 
READ(IOIN,cartlist) 
704 
IF (COORD.NE.'CART') THEN 
705 
Cart=List(1:Nat) 
706 
if (debug) WRITE(*,*) "DBG Path L548 Cart",Cart 
707 
ELSE 
708 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
709 
WRITE(*,*) "I will discard the cartlist" 
710 
Cart=0 
711 
END IF 
712 
END IF 
713 
ELSE 
714 
if (FCart) THEN 
715 
ALLOCATE(Cart(Nat)) 
716 
List=0 
717 
READ(IOIN,cartlist) 
718 
IF (COORD.NE.'CART') THEN 
719 
Cart=List(1:Nat) 
720 
if (debug) WRITE(*,*) "DBG Path L561 Cart",Cart 
721 
ELSE 
722 
WRITE(*,*) "FCart=T is not allowed when Coord='CART'" 
723 
WRITE(*,*) "I will discard the cartlist" 
724 
Cart=0 
725 
END IF 
726 
ELSE 
727 
IF ((PROG=="VASP").OR.(AutoCart)) THEN 
728 
ALLOCATE(Cart(Nat)) 
729 
ELSE 
730 
ALLOCATE(Cart(1)) 
731 
END IF 
732 
Cart=0 
733 
END IF 
734  
735 
if (FFrozen) THEN 
736 
ALLOCATE(Frozen(Nat)) 
737 
List=0 
738 
READ(IOIN,Frozenlist) 
739 
Frozen=List(1:Nat) 
740 
WRITE(*,*) Frozen 
741 
if (debug) WRITE(*,*) "DBG Path L549 Frozen",Frozen 
742 
ELSE 
743 
IF (PROG=="VASP") THEN 
744 
ALLOCATE(Frozen(Nat)) 
745 
ELSE 
746 
ALLOCATE(Frozen(1)) 
747 
END IF 
748 
Frozen=0 
749 
END IF 
750 
END IF 
751  
752  
753 
If (FCart) Call Expand(Cart,NCart,Nat) 
754 
IF (FFrozen) Call Expand(Frozen,NFroz,Nat) 
755  
756 
if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart 
757 
if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen 
758  
759 
if (NMaxPtPath.EQ.1) then 
760 
NMaxPtPath=min(50*NGeomF,2000) 
761 
NMaxPtPath=max(20*NGeomF,NMaxPtPath) 
762 
end if 
763  
764 
! NMaxPtPath=10000 
765  
766 
! We ensure that FTan is in [0.:1.] 
767 
FTan=min(abs(FTan),1.d0) 
768  
769 
! PFL 2011 Mar 14 > 
770 
! Added some printing for analyses with Anapath 
771 
if (debug) THEN 
772 
WRITE(IOOUT,path) 
773 
ELSE 
774 
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC 
775 
! and PATHNAME 
776 
WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat 
777 
WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI 
778 
WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF 
779 
WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC 
780 
WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME) 
781 
END IF 
782  
783 
FTan=1.FTan 
784  
785 
!Thresholds... 
786 
if (SThresh.LE.0) SThresh=0.01d0 
787 
If (GThresh.LE.0) GThresh=SThresh/4. 
788 
SThreshRms=Sthresh*2./3. 
789 
GThreshRms=Gthresh*2./3. 
790  
791 
! First batch of allocations. Arrays that depend only on NGeomI, Nat 
792 
ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) 
793 
ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat)) 
794 
ALLOCATE(MassAt(NAt),OrderInv(Nat)) 
795  
796 
! We read the initial geometries 
797 
Call Read_Geom(input) 
798  
799 
! We convert atome names into atom mass number 
800 
DO I=1,NAt 
801 
! WRITE(*,*) 'I,AtName',I,AtName(I) 
802 
Atome(I)=ConvertNumAt(AtName(I)) 
803 
END DO 
804  
805  
806 
! PFL 23 Sep 2008 > 
807 
! We check that frozen atoms are really frozen, ie that their coordinates do not change 
808 
! between the first geometry and the others. 
809 
IF (NFroz.GT.0) THEN 
810 
ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat)) 
811 
ALLOCATE(x2(nat),y2(nat),z2(nat)) 
812 
ListAtFroz=.FALSE. 
813  
814 
x1(1:Nat)=XyzgeomI(1,1,1:Nat) 
815 
y1(1:Nat)=XyzgeomI(1,2,1:Nat) 
816 
z1(1:Nat)=XyzgeomI(1,3,1:Nat) 
817  
818 
SELECT CASE (NFroz) 
819 
! We have Frozen atoms 
820 
! PFL 28 Oct 2008 > 
821 
! It might happen that the geometries are translated/rotated. 
822 
! So if we have 3 (or more) frozen atoms, we reorient the geometries, based on the first 
823 
! three frozen atoms. 
824  
825 
CASE (3:) 
826 
DO I=1,NFroz 
827 
ListAtFroz(Frozen(I))=.TRUE. 
828 
END DO 
829 
CASE (2) 
830 
IAt1=Frozen(1) 
831 
IAt2=Frozen(2) 
832 
ListAtFroz(Iat1)=.TRUE. 
833 
ListAtFroz(Iat2)=.TRUE. 
834 
x2(1)=x1(Iat2)x1(Iat1) 
835 
y2(1)=y1(Iat2)y1(Iat1) 
836 
z2(1)=z1(Iat2)z1(Iat1) 
837 
Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) 
838 
Dist=1. 
839 
! We scan all atoms to find one that is not aligned with Iat1Iat2 
840 
DO I=1,Nat 
841 
if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle 
842 
x2(2)=x1(Iat2)x1(I) 
843 
y2(2)=y1(Iat2)y1(I) 
844 
z2(2)=z1(Iat2)z1(I) 
845 
Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2) 
846 
ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2) 
847 
if (abs(ca)<0.9) THEN 
848 
IF (Norm2>Dist) THEN 
849 
Iat3=I 
850 
Dist=Norm2 
851 
END IF 
852 
END IF 
853 
END DO 
854 
ListAtFroz(Iat3)=.TRUE. 
855 
if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3 
856 
CASE (1) 
857 
IAt1=Frozen(1) 
858 
ListAtFroz(Iat1)=.TRUE. 
859 
Dist=1. 
860 
! We scan all atoms to find one that is further from At1 
861 
DO I=1,Nat 
862 
if (I.EQ.Iat1) Cycle 
863 
x2(1)=x1(Iat1)x1(I) 
864 
y2(1)=y1(Iat1)y1(I) 
865 
z2(1)=z1(Iat1)z1(I) 
866 
Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) 
867 
IF (Norm1>Dist) THEN 
868 
Iat2=I 
869 
Dist=Norm1 
870 
END IF 
871 
END DO 
872 
IAt2=Frozen(2) 
873 
ListAtFroz(Iat2)=.TRUE. 
874 
if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2 
875 
x2(1)=x1(Iat2)x1(Iat1) 
876 
y2(1)=y1(Iat2)y1(Iat1) 
877 
z2(1)=z1(Iat2)z1(Iat1) 
878 
Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) 
879 
Dist=1. 
880 
! We scan all atoms to find one that is not aligned with Iat1Iat2 
881 
DO I=1,Nat 
882 
if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle 
883 
x2(2)=x1(Iat2)x1(I) 
884 
y2(2)=y1(Iat2)y1(I) 
885 
z2(2)=z1(Iat2)z1(I) 
886 
Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2) 
887 
ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2) 
888 
if (abs(ca)<0.9) THEN 
889 
IF (Norm2>Dist) THEN 
890 
Iat3=I 
891 
Dist=Norm2 
892 
END IF 
893 
END IF 
894 
END DO 
895 
ListAtFroz(Iat3)=.TRUE. 
896 
if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3 
897 
END SELECT 
898  
899 
DO I=2,NGeomI 
900 
! First, we align the Ith geometry on the first one, using only 
901 
! the positions of the "frozen" atoms. (see above: it might be that 
902 
! ListAtFroz contains non frozen atoms usd to align the geometries 
903 
x2(1:Nat)=XyzgeomI(I,1,1:Nat) 
904 
y2(1:Nat)=XyzgeomI(I,2,1:Nat) 
905 
z2(1:Nat)=XyzgeomI(I,3,1:Nat) 
906 
Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot) 
907  
908  
909 
FChkFrozen=.FALSE. 
910 
DO J=1,NFroz 
911 
IAt=Frozen(J) 
912 
IF ((abs(X1(Iat)X2(Iat)).GT.FrozTol).OR. & 
913 
(abs(y1(Iat)y2(Iat)).GT.FrozTol).OR. & 
914 
(abs(z1(Iat)z2(Iat)).GT.FrozTol)) THEN 
915 
IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****" 
916 
FChkFrozen=.TRUE. 
917 
WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change" 
918 
WRITE(*,*) "from geometry 1 to geometry ",I 
919 
END IF 
920 
END DO 
921 
END DO 
922 
IF (FChkFrozen) THEN 
923 
WRITE(*,*) "Some frozen atoms are not frozen... check this and play again" 
924 
STOP 
925 
END IF 
926 
END IF 
927  
928  
929 
N3at=3*Nat 
930 
NFree=N3at6 ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES. 
931  
932 
! if prog=gaussian, we have to read a full input 
933 
! to mimic it later 
934 
if (prog=='GAUSSIAN') THEN 
935 
! We read the Gaussian input file 
936 
! First, the root 
937 
IF (DEBUG) WRITE(*,*) "Reading Gauss Root" 
938 
ALLOCATE(Gauss_Root) 
939 
NULLIFY(Gauss_Root%next) 
940 
Current => Gauss_root 
941 
LineL=1 
942 
DO WHILE (LineL.NE.0) 
943 
READ(IOIN,'(A)') Line 
944 
Line=AdjustL(Line) 
945 
LineL=len(Trim(Line)) 
946 
Idx=INDEX(Line,"chk") 
947 
IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN 
948 
current%Line=TRIM(Line) 
949 
ALLOCATE(current%next) 
950 
Current => Current%next 
951 
Nullify(Current%next) 
952 
END IF 
953 
END DO 
954  
955 
! Current => Gauss_root 
956 
! DO WHILE (ASSOCIATED(Current%next)) 
957 
! WRITE(*,'(1X,A)') Trim(current%line) 
958 
! Current => current%next 
959 
! END DO 
960  
961 
! Now the comment... 
962 
IF (DEBUG) WRITE(*,*) "Reading Gauss Comment" 
963 
ALLOCATE(Gauss_Comment) 
964 
NuLLIFY(Gauss_Comment%Next) 
965 
Current => Gauss_comment 
966 
LineL=1 
967 
DO WHILE (LineL.NE.0) 
968 
READ(IOIN,'(A)') Line 
969 
Line=AdjustL(Line) 
970 
LineL=len(Trim(Line)) 
971 
IF (LineL.NE.0) THEN 
972 
current%Line=TRIM(Line) 
973 
ALLOCATE(current%next) 
974 
Current => Current%next 
975 
Nullify(Current%next) 
976 
END IF 
977 
END DO 
978  
979 
! Current => Gauss_comment 
980 
! DO WHILE (ASSOCIATED(Current%next)) 
981 
! WRITE(*,'(1X,A)') Trim(current%line) 
982 
! Current => current%next 
983 
! END DO 
984  
985 
! Now the charge 
986 
IF (DEBUG) WRITE(*,*) "Reading Gauss Charge" 
987 
READ(IOIN,'(A)') Gauss_Charge 
988 
! WRITE(*,*) Gauss_charge 
989 
! We now read the Paste part... 
990 
ALLOCATE(Gauss_Paste(NAt)) 
991 
LenLine=1 
992 
Iat=0 
993 
Gauss_paste=" " 
994 
DO While (LenLine.GT.0) 
995 
READ(IOIN,'(A)') Line 
996 
Line=AdjustL(Line) 
997 
LenLine=Len_TRIM(Line) 
998 
IF (LenLine.GT.0) Iat=Iat+1 
999 
Idx=Index(Line,'.',BACK=.TRUE.) 
1000 
! WRITE(*,*) Idx,TRIM(Line) 
1001 
! WRITE(*,*) Idx,TRIM(Line(Idx+1:)) 
1002 
Line=ADJUSTL(Line(Idx+1:)) 
1003 
Idx=Index(Line," ") 
1004 
LineL=LEN_TRIM(Line(Idx:)) 
1005 
! WRITE(*,*) LineL,TRIM(Line(Idx:)) 
1006 
IF (LineL.GT.0) THEN 
1007 
Gauss_paste(Iat)=Line(Idx:) 
1008 
END IF 
1009 
END DO 
1010 
IF (Iat.NE.Nat) THEN 
1011 
WRITE(*,*) "I found ", Iat," lines for the geometry instead of ",Nat 
1012 
WRITE(*,*) "ERROR. STOP" 
1013 
WRITE(IOOUT,*) "I found ", Iat," lines for the geometry instead of ",Nat 
1014 
WRITE(IOOUT,*) "ERROR. STOP" 
1015 
STOP 
1016 
END IF 
1017 
! We now read the last part 
1018 
IF (DEBUG) WRITE(*,*) "Reading Gauss End" 
1019 
! READ(IOIN,'(A)') Line 
1020 
ALLOCATE(Gauss_End) 
1021 
NuLLIFY(Gauss_End%Next) 
1022 
Current => Gauss_End 
1023 
LineL=1 
1024 
DO WHILE (1.EQ.1) 
1025 
READ(IOIN,'(A)',END=999) Line 
1026 
Line=AdjustL(Line) 
1027 
LineL=len(Trim(Line)) 
1028 
current%Line=TRIM(Line) 
1029 
ALLOCATE(current%next) 
1030 
Current => Current%next 
1031 
Nullify(Current%next) 
1032 
END DO 
1033 
999 CONTINUE 
1034  
1035 
! ! Write the gaussian input file for testing purposes 
1036 
! Current => Gauss_root 
1037 
! DO WHILE (ASSOCIATED(Current%next)) 
1038 
! WRITE(*,'(1X,A)') Trim(current%line) 
1039 
! Current => current%next 
1040 
! END DO 
1041  
1042 
! WRITE(*,*) 
1043 
! WRITE(*,*) 'Comment original' 
1044  
1045 
! Current => Gauss_comment 
1046 
! DO WHILE (ASSOCIATED(Current%next)) 
1047 
! WRITE(*,'(1X,A)') Trim(current%line) 
1048 
! Current => current%next 
1049 
! END DO 
1050  
1051 
! WRITE(*,*) 
1052 
! WRITE(*,*) Trim(Gauss_charge) 
1053  
1054 
! DO I=1,Nat 
1055 
! WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XyzGeomI(1:3,I,1), TRIM(Gauss_Paste(I)) 
1056 
! END DO 
1057  
1058 
! WRITE(*,*) 
1059 
! Current => Gauss_End 
1060 
! DO WHILE (ASSOCIATED(Current%next)) 
1061 
! WRITE(*,'(1X,A)') Trim(current%line) 
1062 
! Current => current%next 
1063 
! END DO 
1064  
1065 
! WRITE(*,*) 
1066  
1067 
END IF 
1068  
1069 
! if prog=mopac, we have to read a full input to reproduce it later 
1070 
! The structure is: 
1071 
! A MOPAC data set normally consists of one line of keywords, two lines of userdefined text, then the coordinates 
1072 
! Then a blank line or a line of 0. 
1073 
! then the symmetry description. 
1074 
! comment lines start with * and can be anywhere !!! 
1075  
1076 
if (prog=='MOPAC') THEN 
1077 
! We read the Gaussian input file 
1078 
! First, the root 
1079 
IF (DEBUG) WRITE(*,*) "Reading Mopac input" 
1080 
ALLOCATE(Mopac_Root) 
1081 
NULLIFY(Mopac_Root%next) 
1082 
ALLOCATE(Mopac_Comment) 
1083 
NULLIFY(Mopac_Comment%next) 
1084 
ALLOCATE(Mopac_End) 
1085 
NuLLIFY(Mopac_End%Next) 
1086 
Current => Mopac_root 
1087 
CurCom => Mopac_Comment 
1088 
LineL=1 
1089 
NTmp=0 
1090 
DO WHILE (NTmp.LT.3) 
1091 
READ(IOIN,'(A)') Line 
1092 
Line=AdjustL(Line) 
1093 
LineL=len(Trim(Line)) 
1094 
IF (Line(1:1)/="*") THEN 
1095 
IF (NTmp==0) THEN 
1096 
Line2=Line 
1097 
Call UpCase(Line2) 
1098 
Idx=Index(Line2,'GRADIENTS') 
1099 
If (Idx==0) Line=TRIM(Line) // " GRADIENTS" 
1100 
Idx=Index(Line2,'1SCF') 
1101 
If (Idx==0) Line=TRIM(Line) // " 1SCF" 
1102 
END IF 
1103 
current%Line=TRIM(Line) 
1104 
ALLOCATE(current%next) 
1105 
Current => Current%next 
1106 
Nullify(Current%next) 
1107 
NTmp=NTmp+1 
1108 
ELSE 
1109 
CurCom%Line=TRIM(LINE) 
1110 
ALLOCATE(CurCom%Next) 
1111 
CurCom => CurCom%Next 
1112 
NULLIFY(CurCom%Next) 
1113 
END IF 
1114 
END DO 
1115  
1116 
! Current => Mopac_root 
1117 
! DO WHILE (ASSOCIATED(Current%next)) 
1118 
! WRITE(*,'(1X,A)') Trim(current%line) 
1119 
! Current => current%next 
1120 
! END DO 
1121  
1122 
! Now the geometry... that we just skip :) 
1123 
IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry" 
1124 
Mopac_EndGeom="" 
1125 
LineL=1 
1126 
DO WHILE (LineL.NE.0) 
1127 
READ(IOIN,'(A)',END=989) Line 
1128 
Line=AdjustL(Line) 
1129 
LineL=len(Trim(Line)) 
1130 
! The last line might be either blank or filled with 0 
1131 
IF (Line(1:1)=="0") THEN 
1132 
LineL=0 
1133 
Mopac_EndGeom=Trim(Line) 
1134 
END IF 
1135 
IF (Line(1:1)=="*") THEN 
1136 
CurCom%Line=TRIM(LINE) 
1137 
ALLOCATE(CurCom%Next) 
1138 
CurCom => CurCom%Next 
1139 
NULLIFY(CurCom%Next) 
1140 
END IF 
1141 
END DO 
1142  
1143 
! If we are here, there might be something else to read: Mopac_end 
1144  
1145 
! We now read the last part 
1146 
IF (DEBUG) WRITE(*,*) "Reading Mopac End" 
1147 
! READ(IOIN,'(A)') Line 
1148 
Current => Mopac_End 
1149 
LineL=1 
1150 
DO WHILE (1.EQ.1) 
1151 
READ(IOIN,'(A)',END=989) Line 
1152 
Line=AdjustL(Line) 
1153 
LineL=len(Trim(Line)) 
1154 
IF (Line(1:1)/="*") THEN 
1155 
current%Line=TRIM(Line) 
1156 
ALLOCATE(current%next) 
1157 
Current => Current%next 
1158 
Nullify(Current%next) 
1159 
NTmp=NTmp+1 
1160 
ELSE 
1161 
CurCom%Line=TRIM(LINE) 
1162 
ALLOCATE(CurCom%Next) 
1163 
CurCom => CurCom%Next 
1164 
NULLIFY(CurCom%Next) 
1165 
END IF 
1166 
END DO 
1167 
989 CONTINUE 
1168  
1169 
END IF ! IF (PROG="MOPAC") 
1170  
1171 
if ((Prog=="VASP").AND.(input/="VASP")) THEN 
1172 
! Input was not Vasp, so many parameters are missing like lattice 
1173 
! constants... 
1174 
! we read them now ! 
1175 
ALLOCATE(FFF(3,nat)) 
1176 
! First geometry is a bit special for VASP as we have to set 
1177 
! many things 
1178 
IF (DEBUG) WRITE(*,*) "Reading Vasp Parameters" 
1179 
READ(IOIN,'(A)') Vasp_Title 
1180 
READ(IOIN,*) Vasp_param 
1181  
1182 
READ(IOIN,*) Lat_a 
1183 
READ(IOIN,*) Lat_b 
1184 
READ(IOIN,*) Lat_c 
1185  
1186 
Lat_a=Lat_a*Vasp_param 
1187 
Lat_b=Lat_b*Vasp_param 
1188 
Lat_c=Lat_c*Vasp_param 
1189  
1190 
ALLOCATE(NbAtType(nat)) 
1191 
READ(IOIN,'(A)') Vasp_types 
1192 
! First, how many different types ? 
1193 
NbAtType=0 
1194 
READ(Vasp_types,*,END=998) NbAtType 
1195 
998 CONTINUE 
1196 
NbType=0 
1197 
DO WHILE (NbAtType(NbType+1).NE.0) 
1198 
NbType=NbType+1 
1199 
END DO 
1200  
1201 
! Do we know the atom types ? 
1202 
IF (AtTypes(1).EQ.' ') THEN 
1203 
! user has not provided atom types... we have to find them ourselves 
1204 
! by looking into the POTCAR file... 
1205 
INQUIRE(File="POTCAR", EXIST=Tchk) 
1206 
IF (.NOT.Tchk) THEN 
1207 
WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found" 
1208 
STOP 
1209 
END IF 
1210 
OPEN(IOTMP,File="POTCAR") 
1211 
DO I=1,NbType 
1212 
Line='Empty' 
1213 
DO WHILE (Line(1:2).NE.'US') 
1214 
READ(IOTMP,'(A)') Line 
1215 
Line=AdjustL(Line) 
1216 
END DO 
1217 
Line=adjustl(Line(3:)) 
1218 
AtTypes(I)=Line(1:2) 
1219 
END DO 
1220 
if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType) 
1221 
CLOSE(IOTMP) 
1222  
1223 
ELSE !AtTypes(1).EQ.' ' 
1224 
! user has provided atom types 
1225 
NbTypeUser=0 
1226 
DO WHILE (AtTypes(NbTypeUser+1).NE.' ') 
1227 
NbTypeUser=NbTypeUser+1 
1228 
END DO 
1229 
IF (NbType.NE.NbTypeUser) THEN 
1230 
WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes" 
1231 
STOP 
1232 
END IF 
1233 
END IF 
1234  
1235 
IAt=1 
1236 
DO I=1,NbType 
1237 
DO J=1,NbAtType(I) 
1238 
AtName(Iat)=AtTypes(I) 
1239 
Iat=Iat+1 
1240 
END DO 
1241 
END DO 
1242 
DEALLOCATE(NbAtType) 
1243  
1244 
NbTypes=NbType 
1245  
1246 
READ(IOIN,'(A)') Vasp_comment 
1247 
READ(IOIN,'(A)') Vasp_direct 
1248 
V_direct=Adjustl(Vasp_direct) 
1249 
Call UpCase(V_direct) 
1250 
! PFL 2011 Mar 8 > 
1251 
! We have to read the FFF flags : 
1252 
DO I=1,Nat 
1253 
READ(IOIN,*) Xtmp,ytmp,ztmp,FFF(1:3,I) 
1254 
DO J=1,3 
1255 
FFF(J,I)=AdjustL(FFF(J,I)) 
1256 
CALL Upcase(FFF(J,I)) 
1257 
END DO 
1258 
END DO 
1259 
! < PFL 2011 Mar 8 
1260  
1261 
END IF 
1262  
1263 
! In the case of VASP there is always the problem of moving from one side 
1264 
! of the box to the other... 
1265 
! PFL 2011 Mar 14 > 
1266 
! In fact, we have to make this check as soon as either 
1267 
! Input and/or Prog = VASP 
1268 
! PFL 2010 Dec 2 > 
1269 
! Here we deal with input and not prog in fact 
1270 
! if (prog.EQ.'VASP') THEN 
1271 
if ((Input.EQ.'VASP').OR.(Prog.EQ.'VASP')) THEN 
1272 
! < PFL 2010 Dec 2 
1273 
! < PFL 2011 Mar 14 
1274 
Renum=.TRUE. 
1275  
1276 
! V_direct has been set in Read_geom 
1277 
IF (V_direct(1:6).EQ.'DIRECT') THEN 
1278 
Latr(1:3,1)=Lat_a 
1279 
Latr(1:3,2)=Lat_b 
1280 
Latr(1:3,3)=Lat_c 
1281 
B=1. 
1282 
CALL Gaussj(Latr,3,3,B,1,1) 
1283 
ELSE 
1284 
Latr=0. 
1285 
Latr(1,1)=1.d0 
1286 
Latr(2,2)=1.d0 
1287 
Latr(3,3)=1.d0 
1288 
END IF 
1289  
1290 
! Actualization of Frozen using the FFFF... 
1291 
! Frozen has the advantage ie if given, it imposes _ALL_ the FFF flags. 
1292 
IF (Frozen(1).NE.0) THEN 
1293 
WRITE(IOOUT,*) "Frozen is given. Flags of the given POSCAR are overriden" 
1294 
FFF='T' 
1295  
1296 
NFroz=0 
1297 
DO WHILE (Frozen(NFroz+1).NE.0) 
1298 
NFroz=NFroz+1 
1299 
FFF(1:3,Frozen(NFroz))='F' 
1300 
END DO 
1301 
ELSE 
1302 
WRITE(IOOUT,*) "Frozen not given : using Flags of the given POSCAR" 
1303 
NFroz=0 
1304 
Frozen=0 
1305 
DO I=1,Nat 
1306 
IF ((FFF(1,I).EQ.'F').OR.(FFF(2,I).EQ.'F').OR.(FFF(3,I).EQ.'F')) THEN 
1307 
FFF(1:3,I)='F' 
1308 
NFroz=NFroz+1 
1309 
Frozen(NFroz)=I 
1310 
END IF 
1311 
END DO 
1312 
WRITE(IOOUT,*) "Frozen atoms:",Frozen(1:NFroz) 
1313 
END IF 
1314  
1315  
1316 
END IF 
1317  
1318 
IF (Prog=="VASP") THEN 
1319 
IF (Vmd) THEN 
1320 
if (debug) WRITE(*,*) "DBG MAIN, L803, VMD=T,NbTypes,AtTypes",NbTypes,AtTypes(1:NbTypes) 
1321 
Line="" 
1322 
DO I=1,NbTypes 
1323 
Line=TRIM(Line) // " " // TRIM(AdjustL(AtTypes(I))) 
1324 
END DO 
1325 
Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title)) 
1326 
if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title) 
1327 
END IF 
1328 
Call CheckPeriodicBound 
1329 
END IF 
1330  
1331 
IF (COORD=="MIXED") Call TestCart(AutoCart) 
1332  
1333 
SELECT CASE(NFroz) 
1334 
CASE (2) 
1335 
IntFroz=1 
1336 
CASE (3) 
1337 
IntFroz=3 
1338 
CASE (4:) 
1339 
IntFroz=(NFroz2)*3 !(NFroz3)*3+3 
1340 
END SELECT 
1341  
1342 
if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "& 
1343 
,TRIM(COORD),Ncart,NCoord,N3at,NFroz 
1344 
if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart 
1345  
1346 
if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN 
1347 
WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT." 
1348 
COORD="ZMAT" 
1349 
END IF 
1350  
1351 
IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz 
1352  
1353 
SELECT CASE(COORD) 
1354 
CASE ('ZMAT') 
1355 
NCoord=Nfree 
1356 
ALLOCATE(dzdc(3,nat,3,nat)) 
1357 
CASE ('MIXED') 
1358 
SELECT CASE (NCart+NFroz) 
1359 
CASE (1) 
1360 
NCoord=N3At3 
1361 
CASE (2) 
1362 
NCoord=N3At1 
1363 
CASE (3:) 
1364 
NCoord=N3At 
1365 
END SELECT 
1366 
ALLOCATE(dzdc(3,nat,3,nat)) 
1367 
CASE ('BAKER') 
1368 
Symmetry_elimination=0 
1369 
NCoord=3*Nat6Symmetry_elimination !NFree = 3*Nat6 
1370 
ALLOCATE(BMat_BakerT(3*Nat,NCoord)) 
1371 
ALLOCATE(BTransInv(NCoord,3*Nat)) 
1372 
ALLOCATE(BTransInv_local(NCoord,3*Nat)) 
1373 
CASE ('HYBRID') 
1374 
NCoord=N3at 
1375 
CASE ('CART') 
1376 
NCoord=N3at 
1377 
END SELECT 
1378  
1379 
if (debug) WRITE(*,*) "DBG Path, L826, Coord, NCart,NCoord, N3At",TRIM(COORD),Ncart, NCoord,N3at 
1380  
1381 
ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord)) 
1382 
ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat)) 
1383 
ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord)) 
1384 
ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in Path_module.f90 
1385 
ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord)) 
1386 
ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord)) 
1387 
ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in Path_module.f90 
1388 
ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF)) 
1389  
1390 
ZeroVec=0.d0 
1391 
DO I =1, NGeomF 
1392 
IntTangent(I,:)=0.d0 
1393 
END DO 
1394  
1395 
if (debug) THEN 
1396 
Print *, 'L886, IntTangent(1,:)=' 
1397 
Print *, IntTangent(1,:) 
1398 
END IF 
1399  
1400 
if (.NOT.OptReac) Energies(1)=Ereac 
1401 
if (.NOT.OptProd) Energies(NGeomF)=EProd 
1402 
MaxStep=SMax 
1403  
1404 
! We initialize the mass array 
1405 
if (MW) THEN 
1406 
WRITE(*,*) 'Working in MW coordinates' 
1407 
DO I=1,Nat 
1408 
massAt(I)=Mass(Atome(I)) 
1409 
END DO 
1410 
ELSE 
1411 
DO I=1,Nat 
1412 
MassAt(I)=1.d0 
1413 
END DO 
1414 
END IF 
1415  
1416 
WRITE(*,*) "Prog=",TRIM(PROG) 
1417  
1418 
ALLOCATE(FrozAtoms(Nat)) 
1419  
1420 
! IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz 
1421 
! IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen 
1422 
! IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms 
1423 
IF (NFroz.EQ.0) THEN 
1424 
FrozAtoms=.TRUE. 
1425 
ELSE 
1426 
I=1 
1427 
NFr=0 
1428 
FrozAtoms=.False. 
1429 
DO WHILE (Frozen(I).NE.0) 
1430 
IF (Frozen(I).LT.0) THEN 
1431 
DO J=Frozen(I1),abs(Frozen(I)) 
1432 
IF (.NOT.FrozAtoms(J)) THEN 
1433 
NFr=NFr+1 
1434 
FrozAtoms(J)=.TRUE. 
1435 
END IF 
1436 
END DO 
1437 
ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN 
1438 
FrozAtoms(Frozen(I))=.TRUE. 
1439 
NFr=NFr+1 
1440 
END IF 
1441 
I=I+1 
1442 
END DO 
1443 
IF (NFr.NE.NFroz) THEN 
1444 
WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :( ",NFroz,NFr 
1445 
STOP 
1446 
END IF 
1447 
END IF 
1448  
1449 
if (debug) THEN 
1450 
!Print *, 'L968, IntTangent(1,:)=' 
1451 
!Print *, IntTangent(1,:) 
1452 
END IF 
1453  
1454 
! We have read everything we needed... time to work :) 
1455 
IOpt=0 
1456 
FirstTimePathCreate = .TRUE. ! Initialized. 
1457 
Call PathCreate(IOpt) ! IntTangent is also computed in PathCreate. 
1458 
FirstTimePathCreate = .FALSE. 
1459  
1460 
if (debug) THEN 
1461 
Print *, 'L980, IntTangent(1,:)=' 
1462 
Print *, IntTangent(1,:) 
1463 
END IF 
1464  
1465 
! PFL 30 Mar 2008 > 
1466 
! The following is now allocated in Calc_Baker. This is more logical to me 
1467 
! IF (COORD .EQ. "BAKER") Then 
1468 
! ALLOCATE(XprimitiveF(NgeomF,NbCoords)) 
1469 
! ALLOCATE(UMatF(NGeomF,NbCoords,NCoord)) 
1470 
! ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat)) 
1471 
! END IF 
1472 
! < PFL 30 mar 2008 
1473  
1474 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1475 
! 
1476 
! For Debugging purposes 
1477 
! 
1478 
! OptGeom is the index of the geometry which is to be optimized. 
1479 
IF (Optgeom.GT.0) THEN 
1480 
Flag_Opt_Geom=.TRUE. 
1481 
SELECT CASE(Coord) 
1482 
CASE ('CART','HYBRID') 
1483 
!!! CAUTION : PFL 29.AUG.2008 > 
1484 
! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3) 
1485 
! This should be made consistent !!!!!!!!!!!!!!! 
1486 
GeomTmp=Reshape(reshape(XyzGeomF(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
1487 
! GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1488 
!!! < CAUTION : PFL 29.AUG.2008 
1489 
CASE ('ZMAT','MIXED') 
1490 
!GeomTmp=IntCoordF(OptGeom,:) 
1491 
GeomTmp=IntCoordI(OptGeom,:) 
1492 
CASE ('BAKER') 
1493 
!GeomTmp=IntCoordF(OptGeom,:) 
1494 
GeomTmp=IntCoordI(OptGeom,:) 
1495 
CASE DEFAULT 
1496 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP" 
1497 
STOP 
1498 
END SELECT 
1499  
1500 
! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker 
1501 
Flag_Opt_Geom = .TRUE. 
1502 
Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method) 
1503  
1504 
WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom 
1505 
WRITE(*,*) "Before PrintGeom, L1059, Path.f90" 
1506 
Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat) 
1507 
STOP 
1508 
END IF ! IF (Optgeom.GT.0) THEN 
1509  
1510 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1511 
! End of GeomOpt 
1512 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1513  
1514 
IF (PathOnly) THEN 
1515 
WRITE(*,*) "PathOnly=.T. , Stop here" 
1516 
Call Write_path(1) 
1517 
if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar) 
1518 
STOP 
1519 
END IF 
1520  
1521 
if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) 
1522  
1523 
DEALLOCATE(XyzGeomI,IntCoordI) 
1524 
NGeomI=NGeomF 
1525 
ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord)) 
1526  
1527 
IF (DEBUG) WRITE(*,*) ' Back from PathCreate, L965, Path.f90' 
1528 
! we print this path, which is the first one :) called Iteration0 
1529 
! we now have to calculate the gradient for each point (and the energy 
1530 
! if we work at 0K) 
1531  
1532 
if (DEBUG.AND.(COORD.EQ."BAKER")) THEN 
1533 
DO IGeom=1,NGeomF 
1534 
WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom 
1535 
DO I=1,3*Nat6 
1536 
WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) 
1537 
END DO 
1538 
END DO 
1539 
END IF 
1540  
1541 
! To calculate B^T^1 for all extrapolated final geometries: 
1542 
! PFL 18 Jan 2008 > 
1543 
! Added a test for COORD=BAKER before the call to calc_baker_allgeomF 
1544 
IF (COORD.EQ."BAKER") Call Calc_baker_allGeomF() 
1545 
! < PFL 18 Jan 2008 
1546  
1547 
if (DEBUG.AND.(COORD.EQ."BAKER")) THEN 
1548 
DO IGeom=1,NGeomF 
1549 
WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom 
1550 
DO I=1,3*Nat6 
1551 
WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) 
1552 
END DO 
1553 
END DO 
1554 
END IF 
1555  
1556 
IOpt=0 
1557 
Call EgradPath(IOpt,Flag_Opt_Geom) 
1558  
1559 
if (debug) WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt 
1560 
Call Write_path(IOpt) 
1561 
if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) 
1562  
1563 
! We have the geometries, the first gradients... we will generate the first hessian matrices :) 
1564 
! ... v3.71 Only if IniHUp=.TRUE. which is not the default. In this version, we take the hessian 
1565 
! of end points as Unity matrices, and we update them along the path. 
1566  
1567 
! By default, Hess=Id 
1568 
Hess=0. 
1569 
IF(Coord .EQ. "ZMAT") Then 
1570 
! We use the fact that we know that approximate good hessian values 
1571 
! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles 
1572 
DO IGeom=1,NGeomF 
1573 
Hess(IGeom,1,1)=0.5d0 
1574 
Hess(IGeom,2,2)=0.5d0 
1575 
Hess(IGeom,3,3)=0.1d0 
1576 
DO J=1,Nat3 
1577 
Hess(IGeom,3*J+1,3*J+1)=0.5d0 
1578 
Hess(IGeom,3*J+2,3*J+2)=0.1d0 
1579 
Hess(IGeom,3*J+3,3*J+3)=0.02d0 
1580 
END DO 
1581 
IF (HInv) THEN 
1582 
DO I=1,NCoord 
1583 
Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I) 
1584 
END DO 
1585 
END IF 
1586 
END DO 
1587 
ELSE 
1588 
DO I=1,NCoord 
1589 
DO J=1,NGeomF 
1590 
Hess(J,I,I)=1.d0 
1591 
END DO 
1592 
END DO 
1593 
END IF 
1594  
1595 
IF (COORD .EQ. "BAKER") THEN 
1596 
! UMat(NPrim,NCoord) 
1597 
ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat6)) 
1598 
Hprim=0.d0 
1599 
ScanCoord=>Coordinate 
1600 
I=0 
1601 
DO WHILE (Associated(ScanCoord%next)) 
1602 
I=I+1 
1603 
SELECT CASE (ScanCoord%Type) 
1604 
CASE ('BOND') 
1605 
!DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF?? 
1606 
Hprim(I,I) = 0.5d0 
1607 
!END DO 
1608 
CASE ('ANGLE') 
1609 
Hprim(I,I) = 0.2d0 
1610 
CASE ('DIHEDRAL') 
1611 
Hprim(I,I) = 0.1d0 
1612 
END SELECT 
1613 
ScanCoord => ScanCoord%next 
1614 
END DO 
1615  
1616 
! HprimU = H^prim U: 
1617 
HprimU=0.d0 
1618 
DO I=1, 3*Nat6 
1619 
DO J=1,NPrim 
1620 
HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I) 
1621 
END DO 
1622 
END DO 
1623  
1624 
! Hess = U^T HprimU = U^T H^prim U: 
1625 
Hess=0.d0 
1626 
DO I=1, 3*Nat6 
1627 
DO J=1,NPrim 
1628 
! UMat^T is needed below. 
1629 
Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I) 
1630 
END DO 
1631 
END DO 
1632 
DO K=2,NGeomF 
1633 
Hess(K,:,:)=Hess(1,:,:) 
1634 
END DO 
1635 
DEALLOCATE(Hprim,HprimU) 
1636 
end if ! ends IF COORD=BAKER 
1637 
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord)) 
1638 
ALLOCATE(HessTmp(NCoord,NCoord)) 
1639  
1640 
IF (IniHUp) THEN 
1641 
IF (FCalcHess) THEN 
1642 
! The numerical calculation of the Hessian has been switched on 
1643 
DO IGeom=1,NGeomF 
1644 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1645 
GeomTmp=IntCoordF(IGeom,:) 
1646 
ELSE 
1647 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
1648 
END IF 
1649 
Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt) 
1650 
END DO 
1651 
ELSE 
1652 
! We generate a forward hessian and a backward one and we mix them. 
1653 
! First, the forward way: 
1654 
DO IGeom=2,NGeomF1 
1655 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1656 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom1,:) 
1657 
ELSE 
1658 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom1,:,:),(/3*Nat/)) 
1659 
END IF 
1660 
GradTmp=Grad(IGeom1,:)Grad(IGeom,:) 
1661 
HessTmp=Hess(IGeom1,:,:) 
1662 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
1663 
Hess(IGeom,:,:)=HessTmp 
1664 
END DO 
1665  
1666 
! Now, the backward way: 
1667 
HessTmp=Hess(NGeomF,:,:) 
1668 
DO IGeom=NGeomF1,2,1 
1669 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1670 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom+1,:) 
1671 
ELSE 
1672 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) 
1673 
END IF 
1674 
GradTmp=Grad(IGeom,:)Grad(IGeom+1,:) 
1675 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
1676 
alpha=0.5d0*Pi*(IGeom1.)/(NGeomF1.) 
1677 
ca=cos(alpha) 
1678 
sa=sin(alpha) 
1679 
Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:) 
1680 
END DO 
1681 
END IF ! matches IF (FCalcHess) 
1682 
END IF ! matches IF (IniHUp) THEN 
1683  
1684 
! For debug purposes, we diagonalize the Hessian matrices 
1685 
IF (debug) THEN 
1686 
!WRITE(*,*) "DBG Main, L1104,  EigenSystem for Hessian matrices" 
1687 
DO I=1,NGeomF 
1688 
WRITE(*,*) "DBG Main  Point ",I 
1689 
Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord) 
1690 
DO J=1,NCoord 
1691 
WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord)) 
1692 
END DO 
1693 
END DO 
1694 
END IF ! matches IF (debug) THEN 
1695  
1696 
! we have the hessian matrices and the gradients, we can start the optimizations ! 
1697 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1698 
! 
1699 
! Main LOOP : Optimization of the Path 
1700 
! 
1701 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1702 
if (debug) Print *, 'NGeomF=', NGeomF 
1703 
allocation_flag = .TRUE. 
1704  
1705 
Fini=.FALSE. 
1706  
1707 
DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini)) 
1708 
if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt 
1709 
IOpt=IOpt+1 
1710  
1711 
SELECT CASE (COORD) 
1712 
CASE ('ZMAT','MIXED','BAKER') 
1713 
GeomOld_all=IntCoordF 
1714 
CASE ('CART','HYBRID') 
1715 
GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/)) 
1716 
CASE DEFAULT 
1717 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP' 
1718 
STOP 
1719 
END SELECT 
1720  
1721 
IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. & 
1722 
valid("printtangent")) THEN 
1723 
WRITE(*,*) "Visualization of tangents" 
1724 
Idx=min(12,NCoord) 
1725 
DO I=1,NGeomF 
1726 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)0.5d0*IntTangent(I,1:Idx) 
1727 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx) 
1728 
WRITE(*,*) 
1729 
END DO 
1730 
WRITE(*,*) "END of tangents" 
1731 
END IF 
1732  
1733 
IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. & 
1734 
valid("printgrad")) THEN 
1735 
WRITE(*,*) "Visualization of gradients" 
1736 
Idx=min(12,NCoord) 
1737 
DO I=1,NGeomF 
1738 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)1.d0*Grad(I,1:Idx) 
1739 
WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx) 
1740 
WRITe(*,*) 
1741 
END DO 
1742 
WRITE(*,*) "END of gradients" 
1743 
END IF 
1744  
1745 
Fini=.TRUE. 
1746 
IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE. 
1747 
IGeom=1 
1748 
if (debug) WRITE(*,*) '**************************************' 
1749 
if (debug) WRITE(*,*) 'Optimizing reactant' 
1750 
if (debug) WRITE(*,*) '**************************************' 
1751 
SELECT CASE (COORD) 
1752 
CASE ('ZMAT','MIXED','BAKER') 
1753 
SELECT CASE (Step_method) 
1754 
CASE ('RFO') 
1755 
!GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT). 
1756 
Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec) 
1757 
Print *, GradTmp 
1758 
CASE ('GDIIS') 
1759 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).: 
1760 
Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,& 
1761 
Hess(1,:,:),NCoord,allocation_flag,ZeroVec) 
1762 
Print *, 'OptReac, Steps from GDIIS, GeomNew  IntCoordF(1,:)=' 
1763 
Print *, GradTmp 
1764 
CASE ('GEDIIS') 
1765 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).: 
1766 
! Energies are updated in EgradPath.f90 
1767 
Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),& 
1768 
NCoord,allocation_flag,ZeroVec) 
1769 
!Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag) 
1770 
!GradTmp = GeomTmp  IntCoordF(1,:) 
1771 
!Print *, 'Old Geometry:' 
1772 
!Print *, IntCoordF(1,:) 
1773 
Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp  IntCoordF(1,:)=' 
1774 
Print *, GradTmp 
1775 
!Print *, 'GeomTmp:' 
1776 
!Print *, GeomTmp 
1777 
GeomTmp(:) = GradTmp(:) + IntCoordF(1,:) 
1778 
CASE DEFAULT 
1779 
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop" 
1780 
STOP 
1781 
END SELECT 
1782 
CASE ('HYBRID','CART') 
1783 
Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec) 
1784 
CASE DEFAULT 
1785 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP' 
1786 
STOP 
1787 
END SELECT 
1788  
1789 
IF (debug) THEN 
1790 
WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1' 
1791 
Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
1792 
END IF 
1793  
1794 
! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax 
1795 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 
1796 
FactStep=min(1.d0,MaxStep(Igeom)/NormStep) 
1797 
Fini=(Fini.AND.(NormStep.LE.SThresh)) 
1798 
OptReac=(NormStep.GT.SThresh) 
1799 
IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac 
1800 
NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/)))) 
1801 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
1802 
OptReac=(OptReac.OR.(NormGrad.GT.GThresh)) 
1803 
!Print *, 'Grad(1,:):' 
1804 
!Print *, Grad(1,:) 
1805 
IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac 
1806  
1807 
GradTmp=GradTmp*FactStep 
1808  
1809 
! we take care that frozen atoms do not move. 
1810 
IF (NFroz.NE.0) THEN 
1811 
SELECT CASE (COORD) 
1812 
CASE ('ZMAT','MIXED') 
1813 
GradTmp(1:IntFroz)=0.d0 
1814 
CASE ('CART','HYBRID') 
1815 
DO I=1,Nat 
1816 
IF (FrozAtoms(I)) THEN 
1817 
Iat=Order(I) 
1818 
GradTmp(3*Iat2:3*Iat)=0.d0 
1819 
END IF 
1820 
END DO 
1821 
CASE('BAKER') 
1822 
GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..) 
1823 
CASE DEFAULT 
1824 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP' 
1825 
STOP 
1826 
END SELECT 
1827 
END IF ! matches IF (NFroz.NE.0) THEN 
1828  
1829 
IF (debug) THEN 
1830 
!WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1' 
1831 
!Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
1832 
END IF 
1833  
1834 
IF (debug) THEN 
1835 
!WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1' 
1836 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1837 
WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree)) 
1838 
ELSE 
1839 
DO Iat=1,Nat 
1840 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1841 
XyzGeomF(IGeom,1:3,Iat) 
1842 
END DO 
1843 
END IF 
1844 
END IF ! matches IF (debug) THEN 
1845  
1846 
SELECT CASE (COORD) 
1847 
CASE ('ZMAT','MIXED','BAKER') 
1848 
IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1 
1849 
CASE ('HYBRID','CART') 
1850 
XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
1851 
CASE DEFAULT 
1852 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP' 
1853 
STOP 
1854 
END SELECT 
1855  
1856 
IF (debug) THEN 
1857 
WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1' 
1858 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1859 
WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree)) 
1860 
ELSE 
1861 
DO Iat=1,Nat 
1862  
1863 
!!!!!!!!!! There was an OrderInv here... should I put it back ? 
1864 
! It was with indzmat. IndZmat cannot appear here... 
1865 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1866 
XyzGeomF(IGeom,1:3,Iat) 
1867 
END DO 
1868 
END IF 
1869 
END IF ! matches IF (debug) THEN 
1870  
1871 
IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized" 
1872 
ELSE ! Optreac 
1873 
IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini 
1874 
END IF ! matches IF (OptReac) THEN 
1875  
1876 
IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE. 
1877 
IGeom=NGeomF 
1878 
if (debug) WRITE(*,*) '******************' 
1879 
if (debug) WRITE(*,*) 'Optimizing product' 
1880 
if (debug) WRITE(*,*) '******************' 
1881 
SELECT CASE (COORD) 
1882 
CASE ('ZMAT','MIXED','BAKER') 
1883 
Print *, 'Optimizing product' 
1884 
SELECT CASE (Step_method) 
1885 
CASE ('RFO') 
1886 
!GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT). 
1887 
Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec) 
1888 
Print *, GradTmp 
1889 
CASE ('GDIIS') 
1890 
! GradTmp becomes "step" below and has INTENT(OUT): 
1891 
Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),& 
1892 
HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec) 
1893 
Print *, 'OptProd, Steps from GDIIS, GeomNew  IntCoordF(NGeomF,:)=' 
1894 
Print *, GradTmp 
1895 
CASE ('GEDIIS') 
1896 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT): 
1897 
Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),& 
1898 
Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec) 
1899 
Print *, 'OptProd, Steps from GEDIIS, GeomNew  IntCoordF(NGeomF,:)=' 
1900 
Print *, GradTmp 
1901 
CASE DEFAULT 
1902 
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop" 
1903 
STOP 
1904 
END SELECT 
1905 
CASE ('HYBRID','CART') 
1906 
Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec) 
1907 
CASE DEFAULT 
1908 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP' 
1909 
STOP 
1910 
END SELECT 
1911  
1912 
! we have the Newton+RFO step, we will check that the displacement is not greater than Smax 
1913 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 
1914 
FactStep=min(1.d0,MaxStep(IGeom)/NormStep) 
1915 
Fini=(Fini.AND.(NormStep.LE.SThresh)) 
1916 
OptProd=.NOT.(NormStep.LE.SThresh) 
1917 
NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/)))) 
1918 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
1919 
OptProd=(OptProd.OR.(NormGrad.GT.GThresh)) 
1920 
IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd 
1921  
1922 
GradTmp=GradTmp*FactStep 
1923  
1924 
! we take care that frozen atoms do not move 
1925 
IF (NFroz.NE.0) THEN 
1926 
SELECT CASE (COORD) 
1927 
CASE ('ZMAT','MIXED','BAKER') 
1928 
GradTmp(1:IntFroz)=0.d0 
1929 
CASE ('CART','HYBRID') 
1930 
DO I=1,Nat 
1931 
IF (FrozAtoms(I)) THEN 
1932 
Iat=Order(I) 
1933 
GradTmp(3*Iat2:3*Iat)=0.d0 
1934 
END IF 
1935 
END DO 
1936 
CASE DEFAULT 
1937 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP' 
1938 
STOP 
1939 
END SELECT 
1940 
END IF ! matches IF (NFroz.NE.0) THEN 
1941  
1942 
IF (debug) THEN 
1943 
WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt 
1944 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1945 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1946 
ELSE 
1947 
DO Iat=1,Nat 
1948 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1949 
XyzGeomF(IGeom,1:3,Iat) 
1950 
END DO 
1951 
END IF 
1952 
END IF 
1953  
1954 
SELECT CASE (COORD) 
1955 
CASE ('ZMAT','MIXED','BAKER') 
1956 
IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF. 
1957 
CASE ('HYBRID','CART') 
1958 
XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
1959 
CASE DEFAULT 
1960 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP' 
1961 
STOP 
1962 
END SELECT 
1963  
1964 
IF (debug) THEN 
1965 
WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt 
1966 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
1967 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
1968 
ELSE 
1969 
DO Iat=1,Nat 
1970 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
1971 
XyzGeomF(IGeom,1:3,Iat) 
1972 
END DO 
1973 
END IF 
1974 
END IF 
1975 
ELSE ! Optprod 
1976 
IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini 
1977 
END IF !matches IF (OptProd) THEN 
1978  
1979 
IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini 
1980  
1981 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1982 
! 
1983 
! Optimizing other geometries 
1984 
! 
1985 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1986  
1987  
1988  
1989 
DO IGeom=2,NGeomF1 ! matches at L1556 
1990 
if (debug) WRITE(*,*) '****************************' 
1991 
if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom 
1992 
if (debug) WRITE(*,*) '****************************' 
1993  
1994 
SELECT CASE (COORD) 
1995 
CASE ('ZMAT','BAKER','MIXED') 
1996 
GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/)) 
1997 
! PFL 6 Apr 2008 > 
1998 
! Special case, if FTan=0. we do not consider Tangent at all. 
1999 
! To DO: add the full treatment of FTan 
2000 
if (FTan==0.) GradTmp2=ZeroVec 
2001 
! < PFL 6 Apr 2008 
2002 
if (debug) THEN 
2003 
Print *, 'L1500, IntTangent(',IGeom,',:)=' 
2004 
Print *, IntTangent(IGeom,:) 
2005 
END IF 
2006 
!GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT). 
2007 
SELECT CASE (Step_method) 
2008 
CASE ('RFO') 
2009 
Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2) 
2010 
CASE ('GDIIS') 
2011 
!The following has no effect at IOpt==1 
2012 
!Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom 
2013 
Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,& 
2014 
Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2) 
2015 
Print *, 'All others, Steps from GDIIS, GeomNew  IntCoordF(IGeom,:)=' 
2016 
Print *, GradTmp 
2017 
CASE ('GEDIIS') 
2018 
! GradTmp(NCoord) becomes "step" below and has INTENT(OUT): 
2019 
Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),& 
2020 
Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2) 
2021 
Print *, 'All others, Steps from GEDIIS, GeomNew  IntCoordF(IGeom,:)=' 
2022 
Print *, GradTmp 
2023 
CASE DEFAULT 
2024 
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop" 
2025 
STOP 
2026 
END SELECT 
2027 
CASE ('HYBRID','CART') 
2028 
! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat 
2029 
! so we change it 
2030 
DO I=1,Nat 
2031 
DO J=1,3 
2032 
GradTmp2(J+(I1)*3)=XyzTangent(IGeom,(J1)*Nat+I) 
2033 
END DO 
2034 
END DO 
2035 
! GradTmp2=XyzTangent(IGeom,:) 
2036 
! PFL 6 Apr 2008 > 
2037 
! Special case, if FTan=0. we do not consider Tangent at all. 
2038 
! To DO: add the full treatment of FTan 
2039 
if (FTan==0.) GradTmp2=ZeroVec 
2040 
! < PFL 6 Apr 2008 
2041  
2042 
Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2) 
2043 
CASE DEFAULT 
2044 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP' 
2045 
STOP 
2046 
END SELECT 
2047  
2048 
! we take care that frozen atoms do not move 
2049 
IF (NFroz.NE.0) THEN 
2050 
SELECT CASE (COORD) 
2051 
CASE ('ZMAT','MIXED','BAKER') 
2052 
IF (debug) THEN 
2053 
WRITE(*,*) 'Step computed. Coord=',Coord 
2054 
WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz 
2055 
END IF 
2056 
GradTmp(1:IntFroz)=0.d0 
2057 
CASE ('CART','HYBRID') 
2058 
if (debug) THEN 
2059 
WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing" 
2060 
DO I=1,Nat 
2061 
WRITe(*,*) I,GradTmp(3*I2:3*I) 
2062 
END DO 
2063 
END IF 
2064 
DO I=1,Nat 
2065 
IF (FrozAtoms(I)) THEN 
2066 
if (debug) THEN 
2067 
write(*,*) 'Step Computed. Coord=',Coord 
2068 
WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I) 
2069 
END IF 
2070 
Iat=Order(I) 
2071 
GradTmp(3*Iat2:3*Iat)=0.d0 
2072 
END IF 
2073 
END DO 
2074  
2075 
if (debug) THEN 
2076 
WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing" 
2077 
DO I=1,Nat 
2078 
WRITe(*,*) I,GradTmp(3*I2:3*I) 
2079 
END DO 
2080 
END IF 
2081  
2082 
CASE DEFAULT 
2083 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP' 
2084 
STOP 
2085 
END SELECT 
2086 
END IF !matches IF (NFroz.NE.0) THEN 
2087  
2088 
IF (debug) THEN 
2089 
SELECT CASE (COORD) 
2090 
CASE ('ZMAT','MIXED','BAKER') 
2091 
WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree)) 
2092 
WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree)) 
2093 
CASE ('CART','HYBRID') 
2094 
WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree)) 
2095 
WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree)) 
2096 
CASE DEFAULT 
2097 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP' 
2098 
STOP 
2099 
END SELECT 
2100 
END IF ! matches if (debug) THEN 
2101  
2102 
! we check that the maximum displacement is not greater than Smax 
2103 
If (debug) WRITE(*,*) "Fini before test:",Fini 
2104 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 
2105 
FactStep=min(1.d0,maxStep(Igeom)/NormStep) 
2106 
Fini=(Fini.AND.(NormStep.LE.SThresh)) 
2107 
IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini 
2108  
2109 
GradTmp=GraDTmp*FactStep 
2110  
2111 
if (debug) WRITE(*,*) "DBG Main, L1475, FactStep=",FactStep 
2112 
if (debug.AND.(COORD.EQ.'ZMAT')) WRITE(*,'(A,12(1X,F10.6))') 'Scaled Step:',GradTmp(1:min(12,NFree)) 
2113  
2114 
IF (debug) THEN 
2115 
WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt 
2116 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2117 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
2118 
ELSE 
2119 
DO Iat=1,Nat 
2120 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
2121 
XyzGeomF(IGeom,1:3,Iat) 
2122 
END DO 
2123 
END IF 
2124 
END IF ! matches if (debug) THEN 
2125  
2126 
SELECT CASE (COORD) 
2127 
CASE ('ZMAT') 
2128 
IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) 
2129 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
2130 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2131 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) 
2132 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
2133 
OrderInv(indzmat(2,2)),IntcoordF(IGeom,1) 
2134 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
2135 
OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi 
2136 
Idx=4 
2137 
DO Iat=4,Nat 
2138 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
2139 
OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & 
2140 
OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & 
2141 
OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi 
2142 
Idx=Idx+3 
2143 
END DO 
2144  
2145 
CASE ('BAKER') 
2146 
IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) 
2147 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
2148 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2149 
WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull." 
2150  
2151 
CASE ('MIXED') 
2152 
IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) 
2153 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
2154 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2155 
DO Iat=1,NCart 
2156 
WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat2:3*Iat) 
2157 
END DO 
2158 
Idx=3*NCart+1 
2159 
IF (NCart.EQ.1) THEN 
2160 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
2161 
OrderInv(indzmat(2,2)),IntcoordF(IGeom,4) 
2162 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
2163 
OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi 
2164  
2165 
Idx=7 
2166 
END IF 
2167 
IF (NCart.EQ.2) THEN 
2168 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
2169 
OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi 
2170 
Idx=9 
2171 
END IF 
2172  
2173 

2174 
DO Iat=max(NCart+1,4),Nat 
2175 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
2176 
OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & 
2177 
OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & 
2178 
OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi 
2179 
Idx=Idx+3 
2180 
END DO 
2181 
if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) 
2182  
2183 
CASE ('HYBRID','CART') 
2184 
XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) 
2185 
CASE DEFAULT 
2186 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP' 
2187 
STOP 
2188 
END SELECT 
2189  
2190 
IF (debug) THEN 
2191 
WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt 
2192 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2193 
WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) 
2194 
ELSE 
2195 
DO Iat=1,Nat 
2196 
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & 
2197 
XyzGeomF(IGeom,1:3,Iat) 
2198 
END DO 
2199 
END IF 
2200 
END IF ! matches if (debug) THEN 
2201  
2202 
! We project out the tangential part of the gradient to know if we are done 
2203 
GradTmp=Grad(IGeom,:) 
2204 
NormGrad=sqrt(dot_product(GradTmp,GradTmp)) 
2205 
if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad 
2206 
SELECT CASE (COORD) 
2207 
CASE ('ZMAT','MIXED','BAKER') 
2208 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
2209 
Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:)) 
2210 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
2211 
Ca=S/(sqrt(Norm)*NormGrad) 
2212 
S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) 
2213 
GradTmp=GradTmpS/Norm*IntTangent(IGeom,:) 
2214 
CASE ('CART','HYBRID') 
2215 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
2216 
Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:)) 
2217 
Ca=S/(sqrt(Norm)*NormGrad) 
2218 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
2219 
S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) 
2220 
GradTmp=GradTmpS/Norm*XyzTangent(IGeom,:) 
2221 
CASE DEFAULT 
2222 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP' 
2223 
STOP 
2224 
END SELECT 
2225 
IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi 
2226 
NormGrad=sqrt(DOT_Product(GradTmp,GradTmp)) 
2227 
Fini=(Fini.AND.(NormGrad.LE.GThresh)) 
2228 
IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini 
2229  
2230 
END DO ! matches DO IGeom=2,NGeomF1 
2231 
! we save the old gradients 
2232 
GradOld=Grad 
2233 
EPathOld=Energies 
2234  
2235 
! Shall we reparameterize the path ? 
2236 
! PFL 2007/Apr/10 > 
2237 
! We call PathCreate in all cases, it will take care of the 
2238 
! reparameterization, as well as calculating the tangents. 
2239 
! IF (MOD(IOpt,IReparam).EQ.0) THEN 
2240 
XyzGeomI=XyzGeomF 
2241 
IntCoordI=IntCoordF 
2242  
2243 
Call PathCreate(IOpt) 
2244 
! END IF 
2245 
! < PFL 2007/Apr/10 
2246  
2247 
if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)') 'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini 
2248 
IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN 
2249  
2250 
! We have the new path, we calculate its energy and gradient 
2251  
2252 
Call EgradPath(IOpt,Flag_Opt_Geom) 
2253 
!IF(IOPT .EQ. 10) Then 
2254 
! Print *, 'Stopping at my checkpoint.' 
2255 
! STOP !This is my temporary checkpoint. 
2256 
!ENDIF 
2257  
2258 
! PFL 31 Mar 2008 > 
2259 
! Taken from v3.99 but we added a flag... 
2260 
! Added in v3.99 : MaxStep is modified according to the evolution of energies 
2261 
! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2262 
! else it is decreased by 0.8 
2263  
2264 
if (DynMaxStep) THEN 
2265 
if (debug) WRITE(*,*) "Dynamically updating the Maximum step" 
2266 
if (OptReac) THEN 
2267 
If (Energies(1)<EPathOld(1)) Then 
2268 
MaxStep(1)=min(1.1*MaxStep(1),2.*SMax) 
2269 
ELSE 
2270 
MaxStep(1)=max(0.8*MaxStep(1),SMax/2.) 
2271 
END IF 
2272 
END IF 
2273  
2274 
if (OptProd) THEN 
2275 
If (Energies(NGeomF)<EPathOld(NGeomF)) Then 
2276 
MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax) 
2277 
ELSE 
2278 
MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.) 
2279 
END IF 
2280 
END IF 
2281  
2282 
DO IGeom=2,NGeomF1 
2283 
If (Energies(IGeom)<EPathOld(IGeom)) Then 
2284 
MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax) 
2285 
ELSE 
2286 
MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.) 
2287 
END IF 
2288 
END DO 
2289 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
2290 
END IF ! test on DynMaxStep 
2291 
if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) 
2292 
! < PFL 31 MAr 2008 
2293 
! Also XyzGeomF is updated in EgradPath for Baker Case. 
2294 
! It should get updated for other cases too. 
2295  
2296 
DO IGeom=1,NGeomF 
2297 
SELECT CASE (COORD) 
2298 
CASE('ZMAT') 
2299 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2300 
GeomTmp=IntCoordF(IGeom,:) 
2301 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) 
2302 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & 
2303 
OrderInv(indzmat(2,2)),Geomtmp(1) 
2304 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & 
2305 
OrderInv(indzmat(3,2)),Geomtmp(2), & 
2306 
OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi 
2307 
Idx=4 
2308 
DO Iat=4,Nat 
2309 
WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & 
2310 
OrderInv(indzmat(iat,2)),Geomtmp(Idx), & 
2311 
OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, & 
2312 
OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi 
2313 
Idx=Idx+3 
2314 
END DO 
2315 
CASE('BAKER') 
2316 
!WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2317 
!GeomTmp=IntCoordF(IGeom,:) 
2318 
CASE('CART','HYBRID','MIXED') 
2319 
WRITE(IOOUT,'(1X,I5)') Nat 
2320 
WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt 
2321 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2322 
DO I=1,Nat 
2323 
Iat=I 
2324 
If (renum) Iat=OrderInv(I) 
2325 
WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat2:3*Iat) 
2326 
END DO 
2327 
CASE DEFAULT 
2328 
WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP' 
2329 
STOP 
2330 
END SELECT 
2331 
END DO ! matches DO IGeom=1,NGeomF 
2332  
2333 
Call Write_path(IOpt) 
2334  
2335 
! We update the Hessian matrices 
2336 
IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord 
2337 
! First using the displacement from the old path 
2338 
IGeom0=2 
2339 
IGeomF=NGeomF1 
2340 
IF (OptReac) IGeom0=1 
2341 
IF (OptProd) IGeomF=NGeomF 
2342  
2343 
IF (FCalcHess) THEN 
2344 
DO IGeom=IGeom0,IGeomF 
2345 
SELECT CASE (COORD) 
2346 
CASE ('ZMAT','MIXED','BAKER') 
2347 
GeomTmp2=IntCoordF(IGeom,:) 
2348 
CASE ('CART','HYBRID') 
2349 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2350 
CASE DEFAULT 
2351 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
2352 
STOP 
2353 
END SELECT 
2354 
Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt) 
2355 
END DO 
2356 
ELSE 
2357 
DO IGeom=IGeom0,IGeomF 
2358 
GeomTmp=GeomOld_all(IGeom,:) 
2359 
SELECT CASE (COORD) 
2360 
CASE ('ZMAT','MIXED','BAKER') 
2361 
GeomTmp2=IntCoordF(IGeom,:) 
2362 
CASE ('CART','HYBRID') 
2363 
GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
2364 
CASE DEFAULT 
2365 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" 
2366 
STOP 
2367 
END SELECT 
2368 
GeomTmp=GeomTmp2GeomTmp 
2369 
GradTmp=Grad(IGeom,:)GradOld(IGeom,:) 
2370 
HessTmp=Hess(IGeom,:,:) 
2371 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2372 
Hess(IGeom,:,:)=HessTmp 
2373 
END DO ! matches DO IGeom=IGeom0,IGeomF 
2374  
2375 
! Second using the neighbour points 
2376 
IF (HupNeighbour) THEN 
2377 
! Only one point for end points. 
2378 
IF (OptReac) THEN 
2379 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2380 
GeomTmp=IntCoordF(1,:)IntCoordF(2,:) 
2381 
ELSE 
2382 
GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))Reshape(XyzGeomF(2,:,:),(/3*Nat/)) 
2383 
END IF 
2384 
GradTmp=Grad(1,:)Grad(2,:) 
2385 
HessTmp=Hess(1,:,:) 
2386 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2387 
Hess(1,:,:)=HessTmp 
2388 
END IF 
2389  
2390 
IF (OptProd) THEN 
2391 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2392 
GeomTmp=IntCoordF(NGeomF,:)IntCoordF(NGeomF1,:) 
2393 
ELSE 
2394 
GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))Reshape(XyzGeomF(NGeomF1,:,:),(/3*Nat/)) 
2395 
END IF 
2396 
GradTmp=Grad(NGeomF,:)Grad(NGeomF1,:) 
2397 
HessTmp=Hess(NGeomF,:,:) 
2398 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2399 
Hess(NGeomF,:,:)=HessTmp 
2400 
END IF 
2401  
2402 
! Two points for the rest of the path 
2403 
DO IGeom=2,NGeomF1 
2404 
IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN 
2405 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom+1,:) 
2406 
ELSE 
2407 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) 
2408 
END IF 
2409 
GradTmp=Grad(IGeom,:)Grad(IGeom+1,:) 
2410 
HessTmp=Hess(IGeom,:,:) 
2411 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2412  
2413 
IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN 
2414 
GeomTmp=IntCoordF(IGeom,:)IntCoordF(IGeom1,:) 
2415 
ELSE 
2416 
GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))Reshape(XyzGeomF(IGeom1,:,:),(/3*Nat/)) 
2417 
END IF 
2418 
GradTmp=Grad(IGeom,:)Grad(IGeom1,:) 
2419  
2420 
Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) 
2421 
Hess(IGeom,:,:)=HessTmp 
2422 
END DO ! matches DO IGeom=2,NGeomF1 
2423 
END IF !matches IF (HupNeighbour) THEN 
2424 
END IF ! matches IF (FCalcHess) 
2425 
END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc) 
2426 
END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) 
2427  
2428 
IF (PROG=="GAUSSIAN") THEN 
2429 
DEALLOCATE(Gauss_paste) 
2430 
DEALLOCATE(Gauss_root,Gauss_comment,Gauss_end,current) 
2431 
END IF 
2432 
DEALLOCATE(XyzGeomF, IntCoordF) 
2433 
DEALLOCATE(XyzGeomI, IntCoordI) 
2434 
DEALLOCATE(XyzTangent,IntTangent) 
2435 
DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt) 
2436 
DEALLOCATE(GeomOld_all) 
2437  
2438 
if (PROG=="GAUSSIAN") THEN 
2439 
! We deallocate the Gauss_XX lists 
2440 
! transverse the list and deallocate each node 
2441 
previous => Gauss_root ! make current point to head of list 
2442 
DO 
2443 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2444 
current => previous%next ! make list point to next node of head 
2445 
DEALLOCATE(previous) ! deallocate current head node 
2446 
previous => current ! make current point to new head 
2447 
END DO 
2448  
2449 
previous => Gauss_comment ! make current point to head of list 
2450 
DO 
2451 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2452 
current => previous%next ! make list point to next node of head 
2453 
DEALLOCATE(previous) ! deallocate current head node 
2454 
previous => current ! make current point to new head 
2455 
END DO 
2456  
2457 
previous => Gauss_end ! make current point to head of list 
2458 
DO 
2459 
IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer 
2460 
current => previous%next ! make list point to next node of head 
2461 
DEALLOCATE(previous) ! deallocate current head node 
2462 
previous => current ! make current point to new head 
2463 
END DO 
2464  
2465 
DEALLOCATE(current) 
2466 
END IF 
2467  
2468 
DEALLOCATE(GeomTmp, Hess, GradTmp) 
2469 
IF (COORD.EQ.'ZMAT') DEALLOCATE(dzdc) 
2470 
IF (COORD.EQ.'BAKER') THEN 
2471 
DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t) 
2472 
DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local) 
2473 
END IF 
2474  
2475 
WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt 
2476  
2477 
END PROGRAM PathOpt 