## root / src / Path.f90 @ 9

History | View | Annotate | Download (89.5 kB)

1 | 1 | pfleura2 | ! This programs generates and optimizes a reaction path |
---|---|---|---|

2 | 1 | pfleura2 | ! between at least two endpoints. |

3 | 1 | pfleura2 | ! |

4 | 1 | pfleura2 | ! It uses NAMELIST for input, see below ! |

5 | 1 | pfleura2 | ! |

6 | 1 | pfleura2 | ! v3.0 |

7 | 1 | pfleura2 | ! First version entirely in Fortran. |

8 | 1 | pfleura2 | ! It uses routines from Cart (path creation and parameterisation) |

9 | 1 | pfleura2 | ! and from IRC06, especially routines for point optimization. |

10 | 1 | pfleura2 | ! |

11 | 1 | pfleura2 | ! TO DO: |

12 | 1 | pfleura2 | ! 1) Possibility to read and/or calculate some hessian structures |

13 | 1 | pfleura2 | ! 2) Use of internal coordinate to estimate the hessian for stable |

14 | 1 | pfleura2 | ! species in a much better way... |

15 | 1 | pfleura2 | ! |

16 | 1 | pfleura2 | !!!!!!!!!!!!!!! |

17 | 1 | pfleura2 | ! v3.1 |

18 | 1 | pfleura2 | ! The Hessian update includes neighbour points. |

19 | 1 | pfleura2 | ! |

20 | 1 | pfleura2 | ! v3.2 |

21 | 1 | pfleura2 | ! We add the option Hinv that uses the BFGS update of the Hessian inverse |

22 | 1 | pfleura2 | ! |

23 | 1 | pfleura2 | ! v3.3 |

24 | 1 | pfleura2 | ! We add the full option for ZMAT coordinates, where all is done using |

25 | 1 | pfleura2 | ! internal coordinates. |

26 | 1 | pfleura2 | ! v3.31 |

27 | 1 | pfleura2 | ! The step size is controlled using the step norm and not the maximum step |

28 | 1 | pfleura2 | ! component. |

29 | 1 | pfleura2 | ! v3.5 |

30 | 1 | pfleura2 | ! Big modifications of Calc_zmat_const_frag in order to be able |

31 | 1 | pfleura2 | ! to treat VASP geometries. |

32 | 1 | pfleura2 | ! For now, only XYZ geometries are read and written. VASP interface |

33 | 1 | pfleura2 | ! has to be written. |

34 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

35 | 1 | pfleura2 | ! v3.6 |

36 | 1 | pfleura2 | ! Some minor modif in Calc_zmat_const_frag to have nicer programming |

37 | 1 | pfleura2 | ! and a new mode added: Coord=MIXED where part of the system |

38 | 1 | pfleura2 | ! is extrapolated in cartesian coordinates and the rest of it in zmat. |

39 | 1 | pfleura2 | ! it might be Useful for VASP, or when lots of atoms are frozen |

40 | 1 | pfleura2 | ! by default, when Coord=MIXED, all frozen atoms are described in cartesian |

41 | 1 | pfleura2 | ! coordinates. |

42 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

43 | 1 | pfleura2 | ! v3.61 |

44 | 1 | pfleura2 | ! We move further and further to VASP program. New variables added to the 'path' namelist: |

45 | 1 | pfleura2 | ! input: what is the program used for the geometries ? |

46 | 1 | pfleura2 | ! Possible values: xyz (or cart) for Xmol format; VASP. Xyz should be used for Gaussian |

47 | 1 | pfleura2 | ! prog: what is the program used for the calculations ? |

48 | 1 | pfleura2 | ! Possible values for now: gaussian, vasp. |

49 | 1 | pfleura2 | ! For now (02/2007), when prog='vasp' is used, only the initial path is created. |

50 | 1 | pfleura2 | ! That means that we have only added input/output subroutines :-) |

51 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

52 | 1 | pfleura2 | ! v3.7 |

53 | 1 | pfleura2 | ! We try to improve the optimization steps. For that, we first modify many routines |

54 | 1 | pfleura2 | ! so that the new step is calculated in the 'free' degrees of freedom, as it is done |

55 | 1 | pfleura2 | ! in IRC07/ADF for example. That will allow us to impose constraints in an easier way. |

56 | 1 | pfleura2 | ! |

57 | 1 | pfleura2 | ! v3.701 |

58 | 1 | pfleura2 | ! Add a new option for the initial extrapolation of the Hessian |

59 | 1 | pfleura2 | ! |

60 | 1 | pfleura2 | !v3.71 |

61 | 1 | pfleura2 | ! The optimization step has been improved... but the vfree part has not yet been included. |

62 | 1 | pfleura2 | ! This has still to be done :-) |

63 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

64 | 1 | pfleura2 | ! v3.8 |

65 | 1 | pfleura2 | ! I have included part of the vfree procedure: for now freemv returns the Id |

66 | 1 | pfleura2 | ! matrix, that is used to project out the tangent vector. |

67 | 1 | pfleura2 | ! This improves the optimization a lot, but I will have to implement |

68 | 1 | pfleura2 | ! a real freemw subroutine to have constraints there ! |

69 | 1 | pfleura2 | ! v3.81 Technical version |

70 | 1 | pfleura2 | ! I have added a new program 'Ext'. When this one is used, Path creates a Ext.in file, |

71 | 1 | pfleura2 | ! calls Ext.exe file to get a Ext.out file that contains the Energy and the gradient. |

72 | 1 | pfleura2 | ! Format: |

73 | 1 | pfleura2 | ! Ext.in : Xmol format. The second line (comment in Xmol) contains the COORD name. |

74 | 1 | pfleura2 | ! Ext.out: Energy on firt line (1X,D25.15), then gradient in CARTESIAN coordinates |

75 | 1 | pfleura2 | ! (3(1X,D25.15)) format. Natoms lines. |

76 | 1 | pfleura2 | ! TO DO: Gradient could be given in COORD format, ie cartesian for COORD=CART, XYZ or HYBRID |

77 | 1 | pfleura2 | ! ZMAT for CORD=ZMAT (in that case, 6 zeros are there at the begining !) |

78 | 1 | pfleura2 | ! For MIXED, it is equivalent to CART :-) |

79 | 1 | pfleura2 | ! v3.811 |

80 | 1 | pfleura2 | ! I have added a prog="TEST" option, where egradient calls egrad_test subroutine. |

81 | 1 | pfleura2 | ! It then calculates the energy and gradient as it wants and it returns the cartesian |

82 | 1 | pfleura2 | ! gradient. The purpose is to have an analytical approximation of the PES of the system |

83 | 1 | pfleura2 | ! under study to gain some time in testing the program. This is not documented :-) |

84 | 1 | pfleura2 | ! |

85 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

86 | 1 | pfleura2 | ! v3.9 |

87 | 1 | pfleura2 | ! In order to be fully interfaced with VASP, I have to change the architecture |

88 | 1 | pfleura2 | ! of the program. In particular, with VASP, one can use a Para8+NEB routine |

89 | 1 | pfleura2 | ! of VASP to calculate energy+forces for all points of the path. |

90 | 1 | pfleura2 | ! v3.91 |

91 | 1 | pfleura2 | ! One small modification: in the cart list, one can freeze a whole sequence by |

92 | 1 | pfleura2 | ! entering a negative number for the last atom of the sequence |

93 | 1 | pfleura2 | ! for example : cart=1 -160 165 -170 will define atoms from 1 to 160 and from |

94 | 1 | pfleura2 | ! 165 to 170 (included) as cartesian atoms. |

95 | 1 | pfleura2 | ! Of course, Same applies to Frozen |

96 | 1 | pfleura2 | ! v3.92 / v3.93 |

97 | 1 | pfleura2 | ! Slight modifications. One noticeable: when using Vasp, the program analyses the |

98 | 1 | pfleura2 | ! displacements to suggest which atoms could be described in CART. |

99 | 1 | pfleura2 | ! v3.94 |

100 | 1 | pfleura2 | ! VaspMode=Para now implemented. One needs to be careful in defining 'ProgExe'. |

101 | 1 | pfleura2 | ! In particular, one must put the mpirun command in ProgExe !!! There is NO test. |

102 | 1 | pfleura2 | ! v3.95 |

103 | 1 | pfleura2 | ! When using internal coordinates (zmat, hybrid or mixes), Path calculates the number |

104 | 1 | pfleura2 | ! of fragments for each geometry. The reference one (ie the one used to compute the internal |

105 | 1 | pfleura2 | ! coordinates) is now the one with the least fragments. |

106 | 5 | pfleura2 | ! user can also choose one geometry by specifying IGeomRef in the Path namelist. |

107 | 1 | pfleura2 | ! v3.96 |

108 | 1 | pfleura2 | ! I have added a test for the geometry step: it does not allow valence angle to |

109 | 1 | pfleura2 | ! be negative or bigger than 180. I hope that this might help geometry optimization |

110 | 1 | pfleura2 | ! when dealing with nearly linear molecules. |

111 | 1 | pfleura2 | ! v3.961 |

112 | 1 | pfleura2 | ! The options cart and frozen are now LOGICAL: |

113 | 1 | pfleura2 | ! Fcart=T indicates that a &cartlist namelist should be read ! |

114 | 1 | pfleura2 | ! Ffrozen=T indicates that a &frozenlist namelist should be read ! |

115 | 1 | pfleura2 | ! TO DO: Autocart=T indicates that the atoms described in cartesian coordiantes should |

116 | 1 | pfleura2 | ! be found automatically |

117 | 1 | pfleura2 | ! v3.97 |

118 | 1 | pfleura2 | ! Baker has been implemented by P. Dayal: we are interpolating Xint for now. |

119 | 1 | pfleura2 | ! TO do: interpolate Umat ? |

120 | 1 | pfleura2 | ! ----- |

121 | 1 | pfleura2 | ! v3.971 |

122 | 1 | pfleura2 | ! AutoCart=T is being implemented. Goal: having a 'black box' version for Vasp users. |

123 | 1 | pfleura2 | ! Vmd: True if user wants to watch the path using VMD. Only used for VASP for now. |

124 | 1 | pfleura2 | !------------ |

125 | 1 | pfleura2 | ! v3.98 |

126 | 1 | pfleura2 | ! Calc_zmat has been replaced by Calc_zmat_frag which is more robust. |

127 | 1 | pfleura2 | ! Still missing: linear molecules and some tests. |

128 | 1 | pfleura2 | ! |

129 | 1 | pfleura2 | ! v3.981 |

130 | 1 | pfleura2 | ! * Added a non document flag in .valid file that switch on the numerical |

131 | 1 | pfleura2 | ! calculation of the hessian instead of the update. |

132 | 1 | pfleura2 | ! * Added HesUpd variable in the namelist. Default value is MS for path |

133 | 1 | pfleura2 | ! optimization because it does not imposes a positive hessian, and BFGS |

134 | 1 | pfleura2 | ! for geometry optimization |

135 | 1 | pfleura2 | ! v 4.0 |

136 | 1 | pfleura2 | ! * This version has been made public within the CARTE project. |

137 | 1 | pfleura2 | ! * Added: |

138 | 1 | pfleura2 | ! - Step_method to choose the optimization method: RFO or Diis |

139 | 1 | pfleura2 | ! - DynMaxStep: if TRUE, the maximum size of a step is updated at each step. |

140 | 1 | pfleura2 | ! If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. |

141 | 1 | pfleura2 | ! It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0] |

142 | 1 | pfleura2 | ! |

143 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!! |

144 | 1 | pfleura2 | ! v4.1 |

145 | 1 | pfleura2 | ! * Para mode has been partly implemented for Gaussian. |

146 | 1 | pfleura2 | ! VaspMode is thus now 'RunMode' and can be Serial or Para |

147 | 1 | pfleura2 | ! Only Gaussian and VASP are supported for Para mode now. |

148 | 1 | pfleura2 | ! v4.12 |

149 | 1 | pfleura2 | ! Some bugs in unconstrained zmat construction are corrected. |

150 | 1 | pfleura2 | ! v4.13 |

151 | 1 | pfleura2 | ! Prakash has implemented the GEDIIS optimization procedure. |

152 | 1 | pfleura2 | ! v4.14 |

153 | 1 | pfleura2 | ! There is a small problem for some inputs in VASP, where the |

154 | 1 | pfleura2 | ! last geometry is not the one given by the user. |

155 | 1 | pfleura2 | ! It seems to come from the fact that some people try to freeze |

156 | 1 | pfleura2 | ! some unfrozen atoms |

157 | 1 | pfleura2 | ! So some tests have been added to check that frozen atoms do not move. |

158 | 1 | pfleura2 | !!! |

159 | 1 | pfleura2 | ! v4.15 |

160 | 1 | pfleura2 | ! There were some bugs when reading and expanding cartlist and frozenlist |

161 | 1 | pfleura2 | ! I have changed the way frozen atoms are compared by adding partial alignment. |

162 | 1 | pfleura2 | ! v4.16 |

163 | 1 | pfleura2 | ! Some bugs were corrected. |

164 | 1 | pfleura2 | ! v4.17 |

165 | 1 | pfleura2 | ! Added MOPAC as a program. |

166 | 1 | pfleura2 | ! v4.175 |

167 | 1 | pfleura2 | ! Added Three more parameters for defining the input and output names for |

168 | 1 | pfleura2 | ! the calculations. |

169 | 1 | pfleura2 | ! CalcName: Prefix for the files used for the energy and gradient calculations |

170 | 1 | pfleura2 | ! ISuffix: Suffix for the input file |

171 | 1 | pfleura2 | ! OSuffix: suffix for the output file. |

172 | 1 | pfleura2 | ! This maters for Gaussian for example. |

173 | 1 | pfleura2 | ! |

174 | 1 | pfleura2 | ! v4.177 (P) 2010 Nov 22 |

175 | 1 | pfleura2 | ! - We add TurboMole as a new external code. |

176 | 1 | pfleura2 | ! - Subtle change for Input: the default is XYZ whatever the code, except |

177 | 1 | pfleura2 | ! for VASP. |

178 | 1 | pfleura2 | ! |

179 | 1 | pfleura2 | ! v4.178 (P) 2010 Nov 28 |

180 | 1 | pfleura2 | ! - We add the possibility to call CHAMFRE reactive force field using |

181 | 1 | pfleura2 | ! prog=TEST2 or CHAMFRE |

182 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

183 | 1 | pfleura2 | ! OpenPath v1.4 2011 May-June |

184 | 1 | pfleura2 | ! - Also added two flags: CalcEReac and CalcEProd. If .TRUE. Path will |

185 | 1 | pfleura2 | ! compute the energy of the reactant |

186 | 1 | pfleura2 | ! (and/or Product) at the begining of the calculations. |

187 | 1 | pfleura2 | ! This is useful if you start from non optimized points. |

188 | 1 | pfleura2 | ! |

189 | 1 | pfleura2 | ! - Also added more methods for Hessian updates. The following methods are |

190 | 1 | pfleura2 | ! now availables: |

191 | 1 | pfleura2 | ! For inverse Hessian: BFGS, DFP, MS. |

192 | 1 | pfleura2 | ! For Hessian: BFGS, DFP, MS, Bofill |

193 | 1 | pfleura2 | ! |

194 | 1 | pfleura2 | ! -Small change: HesUpd is replaced by Hupdate |

195 | 1 | pfleura2 | ! |

196 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

197 | 1 | pfleura2 | ! TO DO |

198 | 1 | pfleura2 | ! - Trying to improve the output of Path. In particular, we now create |

199 | 1 | pfleura2 | ! a .datl file containing the image index, curvilinear coordinate and energy |

200 | 1 | pfleura2 | ! in the spirit of the AnaPath tool. |

201 | 1 | pfleura2 | ! AnaPath is still usefull as it can extract geometrical parameters on top |

202 | 1 | pfleura2 | ! of the energy |

203 | 1 | pfleura2 | ! - When using Hessian update, it is more efficient to use Cholesky |

204 | 1 | pfleura2 | ! decomposition, and to update this decomposition rather than the Hessian. |

205 | 1 | pfleura2 | ! See Nocedal&Wright p141. |

206 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

207 | 5 | pfleura2 | ! |

208 | 1 | pfleura2 | ! Project Name has changed to Opt'n Path |

209 | 1 | pfleura2 | ! OpenPaht 2011 May-June -> Opt'n Path v1.4 |

210 | 5 | pfleura2 | ! |

211 | 5 | pfleura2 | ! 2013 Jan 18 Opt'n Path v1.47 (going to 1.5 slowly) |

212 | 5 | pfleura2 | ! We add Siesta (and start thinking of add CP2K) |

213 | 5 | pfleura2 | ! as a new "energy engine" code |

214 | 5 | pfleura2 | ! TO DO to go to 1.5: finish cleaning... |

215 | 7 | pfleura2 | ! |

216 | 5 | pfleura2 | ! 2013 Feb Opt'n Path v1.48 |

217 | 5 | pfleura2 | ! We add some keyword for more output for Geometry Optimizations |

218 | 5 | pfleura2 | ! GeomFile: name of the file to print out the geometries and their energy |

219 | 5 | pfleura2 | ! as Xyz file. (only format for now) |

220 | 7 | pfleura2 | ! |

221 | 7 | pfleura2 | ! 2013 Feb opt'n Path v1.49 |

222 | 7 | pfleura2 | ! We add the possibility to analyze the geometries in terms |

223 | 7 | pfleura2 | ! of Bonds, Angles, dihedrals |

224 | 7 | pfleura2 | ! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path) |

225 | 7 | pfleura2 | ! This is basicaly the merging of AnaPath into Opt'n Path... at last ! |

226 | 7 | pfleura2 | ! This is done by adding a Flag into the Path namelist: |

227 | 7 | pfleura2 | ! AnaGeom: logical, if true Path looks for the AnaList namelist |

228 | 7 | pfleura2 | ! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path |

229 | 7 | pfleura2 | ! will read it and print the values of the variable in a Path.dat file |

230 | 7 | pfleura2 | ! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. |

231 | 1 | pfleura2 | |

232 | 1 | pfleura2 | PROGRAM PathOpt |

233 | 1 | pfleura2 | |

234 | 1 | pfleura2 | use VarTypes |

235 | 1 | pfleura2 | use Path_module |

236 | 1 | pfleura2 | use Io_module |

237 | 1 | pfleura2 | |

238 | 1 | pfleura2 | IMPLICIT NONE |

239 | 1 | pfleura2 | |

240 | 5 | pfleura2 | CHARACTER(132) :: FileIn, FileOut, Line |

241 | 1 | pfleura2 | |

242 | 1 | pfleura2 | INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF |

243 | 1 | pfleura2 | INTEGER(KINT) :: IOpt |

244 | 5 | pfleura2 | INTEGER(KINT) :: Idx |

245 | 1 | pfleura2 | INTEGER(KINT) :: N3at, NFree, Nfr |

246 | 1 | pfleura2 | |

247 | 1 | pfleura2 | INTEGER(KINT) :: List(2*MaxFroz) |

248 | 1 | pfleura2 | |

249 | 5 | pfleura2 | REAL(KREAL) :: E,FactStep !E=Energy |

250 | 1 | pfleura2 | REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad |

251 | 1 | pfleura2 | REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E |

252 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord |

253 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord) |

254 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord |

255 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: HessTmp(:,:) |

256 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim) |

257 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6)) |

258 | 1 | pfleura2 | LOGICAL, SAVE :: allocation_flag |

259 | 1 | pfleura2 | |

260 | 1 | pfleura2 | ! Flag to see if frozen atoms are frozen (see below) |

261 | 1 | pfleura2 | LOGICAL :: FChkFrozen |

262 | 1 | pfleura2 | |

263 | 1 | pfleura2 | ! Energies for all points along the path at the previous iteration |

264 | 7 | pfleura2 | REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF |

265 | 1 | pfleura2 | ! Maximum step allowed for this geometry |

266 | 7 | pfleura2 | REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF |

267 | 1 | pfleura2 | |

268 | 1 | pfleura2 | ! these are used to read temporary coordinates |

269 | 1 | pfleura2 | LOGICAL :: FFrozen,FCart |

270 | 1 | pfleura2 | |

271 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

272 | 1 | pfleura2 | ! |

273 | 1 | pfleura2 | ! To analyse the frozen atoms |

274 | 1 | pfleura2 | ! |

275 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat |

276 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat |

277 | 1 | pfleura2 | REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist |

278 | 1 | pfleura2 | LOGICAL, ALLOCATABLE :: ListAtFroz(:) |

279 | 1 | pfleura2 | INTEGER(KINT) :: Iat1, Iat2, Iat3 |

280 | 1 | pfleura2 | |

281 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

282 | 1 | pfleura2 | |

283 | 1 | pfleura2 | LOGICAL :: Debug, Fini,Flag_Opt_Geom |

284 | 1 | pfleura2 | |

285 | 1 | pfleura2 | ! INTEGER(KINT), EXTERNAL :: Iargc |

286 | 1 | pfleura2 | INTEGER(KINT), EXTERNAL :: ConvertNumAt |

287 | 7 | pfleura2 | INTEGER(KINT) :: IoS |

288 | 1 | pfleura2 | |

289 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

290 | 1 | pfleura2 | ! |

291 | 1 | pfleura2 | ! For Debugging purposes |

292 | 1 | pfleura2 | ! |

293 | 1 | pfleura2 | LOGICAL :: FCalcHess |

294 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

295 | 1 | pfleura2 | |

296 | 1 | pfleura2 | |

297 | 1 | pfleura2 | INTERFACE |

298 | 1 | pfleura2 | function valid(string) result (isValid) |

299 | 1 | pfleura2 | CHARACTER(*), intent(in) :: string |

300 | 1 | pfleura2 | logical :: isValid |

301 | 1 | pfleura2 | END function VALID |

302 | 1 | pfleura2 | |

303 | 1 | pfleura2 | SUBROUTINE Read_Geom(Input) |

304 | 1 | pfleura2 | CHARACTER(32), INTENT(IN) :: input |

305 | 1 | pfleura2 | END SUBROUTINE READ_GEOM |

306 | 1 | pfleura2 | |

307 | 5 | pfleura2 | |

308 | 1 | pfleura2 | subroutine hupdate_all (n, dx, dgrad, HessOld) |

309 | 1 | pfleura2 | |

310 | 1 | pfleura2 | IMPLICIT NONE |

311 | 1 | pfleura2 | |

312 | 1 | pfleura2 | INTEGER, PARAMETER :: KINT=KIND(1) |

313 | 1 | pfleura2 | INTEGER, PARAMETER :: KREAL=KIND(1.0D0) |

314 | 1 | pfleura2 | |

315 | 1 | pfleura2 | INTEGER(KINT), INTENT(IN) :: n |

316 | 1 | pfleura2 | real(KREAL) :: dx(n), dgrad(n), HessOld(n*n) |

317 | 1 | pfleura2 | END subroutine hupdate_all |

318 | 1 | pfleura2 | |

319 | 1 | pfleura2 | SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess) |

320 | 1 | pfleura2 | |

321 | 1 | pfleura2 | IMPLICIT NONE |

322 | 1 | pfleura2 | |

323 | 1 | pfleura2 | INTEGER, PARAMETER :: KINT=KIND(1) |

324 | 1 | pfleura2 | INTEGER, PARAMETER :: KREAL=KIND(1.D0) |

325 | 1 | pfleura2 | |

326 | 1 | pfleura2 | INTEGER(KINT), INTENT(IN) :: NFREE |

327 | 1 | pfleura2 | REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree) |

328 | 1 | pfleura2 | REAL(KREAL), INTENT(IN) :: ds(Nfree) |

329 | 1 | pfleura2 | REAL(KREAL), INTENT(IN) :: dGrad(Nfree) |

330 | 1 | pfleura2 | |

331 | 1 | pfleura2 | END subroutine Hinvup_BFGS_new |

332 | 1 | pfleura2 | |

333 | 1 | pfleura2 | |

334 | 8 | pfleura2 | FUNCTION InString(Line,String,Case,Clean,Back) Result(Pos) |

335 | 8 | pfleura2 | |

336 | 8 | pfleura2 | Use VarTypes |

337 | 8 | pfleura2 | |

338 | 8 | pfleura2 | implicit none |

339 | 8 | pfleura2 | ! Input |

340 | 8 | pfleura2 | CHARACTER(*), INTENT(IN) :: Line |

341 | 8 | pfleura2 | CHARACTER(*), INTENT(IN) :: String |

342 | 8 | pfleura2 | LOGICAL, OPTIONAL, INTENT(IN) :: Case |

343 | 8 | pfleura2 | LOGICAL, OPTIONAL, INTENT(IN) :: Back |

344 | 8 | pfleura2 | CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean |

345 | 8 | pfleura2 | |

346 | 8 | pfleura2 | ! Output |

347 | 8 | pfleura2 | ! the position of String in Line (the first one) unless Back is present |

348 | 8 | pfleura2 | INTEGER(KINT) :: Pos |

349 | 8 | pfleura2 | END FUNCTION InString |

350 | 8 | pfleura2 | |

351 | 8 | pfleura2 | SUBROUTINE die(routine, msg, file, line, unit) |

352 | 8 | pfleura2 | |

353 | 8 | pfleura2 | Use VarTypes |

354 | 8 | pfleura2 | Use io_module |

355 | 8 | pfleura2 | |

356 | 8 | pfleura2 | implicit none |

357 | 8 | pfleura2 | character(len=*), intent(in) :: routine, msg |

358 | 8 | pfleura2 | character(len=*), intent(in), optional :: file |

359 | 8 | pfleura2 | integer(KINT), intent(in), optional :: line, unit |

360 | 8 | pfleura2 | |

361 | 8 | pfleura2 | END SUBROUTINE die |

362 | 8 | pfleura2 | |

363 | 8 | pfleura2 | SUBROUTINE Warning(routine, msg, file, line, unit) |

364 | 8 | pfleura2 | |

365 | 8 | pfleura2 | Use VarTypes |

366 | 8 | pfleura2 | Use io_module |

367 | 8 | pfleura2 | |

368 | 8 | pfleura2 | implicit none |

369 | 8 | pfleura2 | |

370 | 8 | pfleura2 | character(len=*), intent(in) :: routine, msg |

371 | 8 | pfleura2 | character(len=*), intent(in), optional :: file |

372 | 8 | pfleura2 | integer(KINT), intent(in), optional :: line, unit |

373 | 8 | pfleura2 | |

374 | 8 | pfleura2 | END SUBROUTINE Warning |

375 | 8 | pfleura2 | |

376 | 8 | pfleura2 | |

377 | 1 | pfleura2 | END INTERFACE |

378 | 1 | pfleura2 | |

379 | 1 | pfleura2 | Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, & |

380 | 1 | pfleura2 | MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, & |

381 | 1 | pfleura2 | SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, & |

382 | 1 | pfleura2 | input,prog, ProgExe,RunMode, AtTypes,poscar,PathName, & |

383 | 1 | pfleura2 | OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, & |

384 | 1 | pfleura2 | FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, & |

385 | 1 | pfleura2 | NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, & |

386 | 8 | pfleura2 | NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile |

387 | 1 | pfleura2 | |

388 | 1 | pfleura2 | Namelist/cartlist/list |

389 | 1 | pfleura2 | Namelist/frozenlist/list |

390 | 7 | pfleura2 | Namelist/analist/nb |

391 | 1 | pfleura2 | |

392 | 1 | pfleura2 | |

393 | 1 | pfleura2 | Flag_Opt_Geom = .FALSE. ! Initialized. |

394 | 1 | pfleura2 | |

395 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!! |

396 | 1 | pfleura2 | ! |

397 | 1 | pfleura2 | ! We print the Version info in file |

398 | 1 | pfleura2 | ! |

399 | 5 | pfleura2 | WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013" |

400 | 1 | pfleura2 | |

401 | 8 | pfleura2 | ! We read the files names |

402 | 1 | pfleura2 | ! SELECT CASE (iargc()) |

403 | 1 | pfleura2 | SELECT CASE (command_argument_count()) |

404 | 1 | pfleura2 | CASE (2) |

405 | 1 | pfleura2 | call getarg(1,FileIn) |

406 | 1 | pfleura2 | OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') |

407 | 1 | pfleura2 | call getarg(2,FileOut) |

408 | 1 | pfleura2 | OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN') |

409 | 1 | pfleura2 | CASE (1) |

410 | 1 | pfleura2 | call getarg(1,FileIn) |

411 | 1 | pfleura2 | IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN |

412 | 9 | pfleura2 | WRITE(*,*) "Opt'n Path mini-help" |

413 | 1 | pfleura2 | WRITE(*,*) "--------------" |

414 | 1 | pfleura2 | WRITE(*,*) "" |

415 | 1 | pfleura2 | WRITE(*,*) "Use: Path Input_file Output_file" |

416 | 1 | pfleura2 | WRITE(*,*) "Input_file starts with a Namelist called path" |

417 | 1 | pfleura2 | WRITE(*,*) "" |

418 | 9 | pfleura2 | ! The following has been generated (March 19, 2013) using the Mini_help.tex file |

419 | 9 | pfleura2 | ! 1) Edit Mini_help.tex |

420 | 9 | pfleura2 | ! 2) latex Mini_help.tex |

421 | 9 | pfleura2 | ! 3) catdvi -e 1 -U Mini_help.dvi | sed -re "s/\[U\+2022\]/-/g" | sed -re "s/([^^[:space:]])\s+/\1 /g" > file.txt |

422 | 9 | pfleura2 | ! 4) edit file.txt to keep only the following lines, and to reformat a bit |

423 | 9 | pfleura2 | ! 5) awk '{print "WRITE(*,*) \"" $0 "\""}' file.txt > help |

424 | 9 | pfleura2 | ! 5) cut&paste help in Path.f90 ... |

425 | 9 | pfleura2 | WRITE(*,*) " &path" |

426 | 9 | pfleura2 | WRITE(*,*) " nat=3, ! Number of atoms" |

427 | 9 | pfleura2 | WRITE(*,*) " ngeomi=3, ! Number of initial geometries" |

428 | 9 | pfleura2 | WRITE(*,*) " ngeomf=12, !Number of geometries along the path" |

429 | 9 | pfleura2 | WRITE(*,*) " OptReac=.T., ! Do you want to optimize the reactants ?" |

430 | 9 | pfleura2 | WRITE(*,*) " OptProd=.T., ! Optimize the products" |

431 | 9 | pfleura2 | WRITE(*,*) " coord='zmat', ! We use Z-matrix coordinates" |

432 | 9 | pfleura2 | WRITE(*,*) " maxcyc=31, ! Max number of iterations" |

433 | 9 | pfleura2 | WRITE(*,*) " IReparam=2,! re-distribution of points along the path every 2 iterations" |

434 | 9 | pfleura2 | WRITE(*,*) " ISpline=50, ! Start using spline interpolation at iteration 50" |

435 | 9 | pfleura2 | WRITE(*,*) " Hinv=.T. , ! Use inverse of the Hessian internally (default: T)" |

436 | 9 | pfleura2 | WRITE(*,*) " MW=T, ! Works in Mass Weighted coordiante (default T)" |

437 | 9 | pfleura2 | WRITE(*,*) " PathName='Path_HCN_zmat_test', ! Name of the file used for path outputs" |

438 | 9 | pfleura2 | WRITE(*,*) " prog='gaussian',! we use G03 to get energy and gradients" |

439 | 9 | pfleura2 | WRITE(*,*) " SMax=0.1 ! Displacement cannot exceed 0.1 atomic units (or mass weighted at. unit)" |

440 | 9 | pfleura2 | WRITE(*,*) " /" |

441 | 9 | pfleura2 | WRITE(*,*) " 3" |

442 | 9 | pfleura2 | WRITE(*,*) " Energy : 0.04937364" |

443 | 9 | pfleura2 | WRITE(*,*) " H 0.0000 0.0000 0.0340" |

444 | 9 | pfleura2 | WRITE(*,*) " C 0.0000 0.0000 1.1030" |

445 | 9 | pfleura2 | WRITE(*,*) " N 0.0000 0.0000 2.2631" |

446 | 9 | pfleura2 | WRITE(*,*) " 3" |

447 | 9 | pfleura2 | WRITE(*,*) " Energy : 0.04937364" |

448 | 9 | pfleura2 | WRITE(*,*) " H 0.0000 1.1000 1.1030" |

449 | 9 | pfleura2 | WRITE(*,*) " C 0.0000 0.0000 1.1030" |

450 | 9 | pfleura2 | WRITE(*,*) " N 0.0000 0.0000 2.2631" |

451 | 9 | pfleura2 | WRITE(*,*) "3" |

452 | 9 | pfleura2 | WRITE(*,*) " CNH" |

453 | 9 | pfleura2 | WRITE(*,*) "H 0.000000 0.000000 3.3" |

454 | 9 | pfleura2 | WRITE(*,*) "C 0.000000 0.000000 1.1" |

455 | 9 | pfleura2 | WRITE(*,*) "N 0.000000 0.000000 2.26" |

456 | 9 | pfleura2 | WRITE(*,*) "%chk=Test" |

457 | 9 | pfleura2 | WRITE(*,*) "#P AM1 FORCE" |

458 | 9 | pfleura2 | WRITE(*,*) " HCN est bien" |

459 | 9 | pfleura2 | WRITE(*,*) "" |

460 | 9 | pfleura2 | WRITE(*,*) "0,1" |

461 | 9 | pfleura2 | WRITE(*,*) "H 0.000000 0.000000 0.000000" |

462 | 9 | pfleura2 | WRITE(*,*) "C 0.000000 0.000000 1.000" |

463 | 9 | pfleura2 | WRITE(*,*) "N 0.000000 0.000000 3.00" |

464 | 9 | pfleura2 | WRITE(*,*) "" |

465 | 9 | pfleura2 | WRITE(*,*) "*******" |

466 | 9 | pfleura2 | WRITE(*,*) "* Compulsory variables are:" |

467 | 9 | pfleura2 | WRITE(*,*) "*******" |

468 | 9 | pfleura2 | WRITE(*,*) "NGeomi: Number of geometries defining the Initial path" |

469 | 9 | pfleura2 | WRITE(*,*) "NGeomf: Number of geometries defining the Final path" |

470 | 9 | pfleura2 | WRITE(*,*) "Nat: Number of atoms" |

471 | 9 | pfleura2 | WRITE(*,*) "" |

472 | 9 | pfleura2 | WRITE(*,*) "*******" |

473 | 9 | pfleura2 | WRITE(*,*) "* Other options:" |

474 | 9 | pfleura2 | WRITE(*,*) "*******" |

475 | 9 | pfleura2 | WRITE(*,*) "Input: String that indicates the type of the input geometries. Accepted values are: Cart (or" |

476 | 9 | pfleura2 | WRITE(*,*) " Xmol or Xyz), Vasp, Turbomole or Siesta." |

477 | 9 | pfleura2 | WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations." |

478 | 9 | pfleura2 | WRITE(*,*) " Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext." |

479 | 9 | pfleura2 | WRITE(*,*) "" |

480 | 9 | pfleura2 | WRITE(*,*) "" |

481 | 9 | pfleura2 | WRITE(*,*) " - In case of a Gaussian calculations, input must be set to Cart. One example of a gaus-" |

482 | 9 | pfleura2 | WRITE(*,*) " sian input should be added at the end of the input file.See example file Test_HCN_zmat_g03.path." |

483 | 9 | pfleura2 | WRITE(*,*) " - In the case of a VASP calculation, if input is set to Cart, then the preamble of a" |

484 | 9 | pfleura2 | WRITE(*,*) " VASP calculation s" |

485 | 9 | pfleura2 | WRITE(*,*) " - Mopac: Examples using sequential call to MOPAC2009 to calculate the energies and" |

486 | 9 | pfleura2 | WRITE(*,*) " forces along the path. hould be added at the end of the input file. See example file" |

487 | 9 | pfleura2 | WRITE(*,*) " Test_VASP_cart.path. In the case of a VASP calculation, one should also give value" |

488 | 9 | pfleura2 | WRITE(*,*) " of the RunMode variable ." |

489 | 9 | pfleura2 | WRITE(*,*) " - In the case of a SIESTA calculation, an example of a Siesta input file should be added" |

490 | 9 | pfleura2 | WRITE(*,*) " at the end of the input file. See Test_Siesta.path." |

491 | 9 | pfleura2 | WRITE(*,*) "" |

492 | 9 | pfleura2 | WRITE(*,*) "RunMode: When running on a multi-processor machine, this indicates wether Opt'n Path" |

493 | 9 | pfleura2 | WRITE(*,*) " should calculate the energy and gradient of the whole path in parallel or not. User has" |

494 | 9 | pfleura2 | WRITE(*,*) " two options. One is to calculate the energy and gradient of each point sequentially. This" |

495 | 9 | pfleura2 | WRITE(*,*) " is usefull when running on one (or two) processors. In this case, RunMode should be put" |

496 | 9 | pfleura2 | WRITE(*,*) " to SERIAL.When running in parallel with 8 or more processors, one can use VASP to" |

497 | 9 | pfleura2 | WRITE(*,*) " calculate simultaneously the energies and gradients for all points, as in a normal NEB cal-" |

498 | 9 | pfleura2 | WRITE(*,*) " culation. In this case, RunMode must be set to PARA. For now, this is usefull only for VASP." |

499 | 9 | pfleura2 | WRITE(*,*) "" |

500 | 9 | pfleura2 | WRITE(*,*) "ProgExe: Name (with full path) of the executable to be used to get energies and gradients." |

501 | 9 | pfleura2 | WRITE(*,*) " For example, if VASP is used in parallel, one might have something like:" |

502 | 9 | pfleura2 | WRITE(*,*) " ProgExe='/usr/local/mpich/bin/mpirun -machinefile machine -np 8 ~/bin/VASP_46'." |

503 | 9 | pfleura2 | WRITE(*,*) " Another option that I use, is to put ProgExe='./run_vasp' and I put every option needed" |

504 | 9 | pfleura2 | WRITE(*,*) " to run VASP into the run_vasp file." |

505 | 9 | pfleura2 | WRITE(*,*) "EReac: (REAL) By default, Opt'n Path does not compute the energy of the reactants and" |

506 | 9 | pfleura2 | WRITE(*,*) " products. This thus indicates the reactants energy to Opt'n Path to have better plots at" |

507 | 9 | pfleura2 | WRITE(*,*) " the end." |

508 | 9 | pfleura2 | WRITE(*,*) "EProd: (REAL) By default, Opt'n Path does not compute the energy of the reactants and" |

509 | 9 | pfleura2 | WRITE(*,*) " products. This thus indicates the products energy to Opt'n Path to have better plots." |

510 | 9 | pfleura2 | WRITE(*,*) "" |

511 | 9 | pfleura2 | WRITE(*,*) "PathName: Prefix used to save the path. Default is Path." |

512 | 9 | pfleura2 | WRITE(*,*) "Poscar: string that will be used as the prefix for the different POSCAR files in a VASP calcu-" |

513 | 9 | pfleura2 | WRITE(*,*) " lations. Usefull only if PathOnly=.TRUE., not used for internal calculations." |

514 | 9 | pfleura2 | WRITE(*,*) "CalcName: Prefix for the files used for the energy and gradient calculations" |

515 | 9 | pfleura2 | WRITE(*,*) "" |

516 | 9 | pfleura2 | WRITE(*,*) "ISuffix: Suffix for the input file used for energy and gradient calculations. The full inputfile" |

517 | 9 | pfleura2 | WRITE(*,*) " name will thus be CalcName.ISuffix." |

518 | 9 | pfleura2 | WRITE(*,*) "OSuffix: Suffix for the output file used for energy and gradient calculations. The full outputfile" |

519 | 9 | pfleura2 | WRITE(*,*) " name will thus be CalcName.OSuffix." |

520 | 9 | pfleura2 | WRITE(*,*) "" |

521 | 9 | pfleura2 | WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates. Valid only for" |

522 | 9 | pfleura2 | WRITE(*,*) " Coord=Zmat, Hybrid or Mixed." |

523 | 9 | pfleura2 | WRITE(*,*) "Fact: REAL used to define if two atoms are linked. If d(A,B) =< fact* (rcov(A) + rcov(B))," |

524 | 9 | pfleura2 | WRITE(*,*) " then A and B are considered Linked." |

525 | 9 | pfleura2 | WRITE(*,*) "" |

526 | 9 | pfleura2 | WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info." |

527 | 9 | pfleura2 | WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:" |

528 | 9 | pfleura2 | WRITE(*,*) " - CART (or Xyz): works in cartesian" |

529 | 9 | pfleura2 | WRITE(*,*) " - Zmat: works in internal coordinates (Zmat)" |

530 | 9 | pfleura2 | WRITE(*,*) " - Mixed: frozen atoms, as well as atoms defined by the 'cart' array(see below) are" |

531 | 9 | pfleura2 | WRITE(*,*) " describe in CARTESIAN, whereas the others are described in Zmat" |

532 | 9 | pfleura2 | WRITE(*,*) " - Baker: use of Baker coordinates, also called delocalized internal coordinates" |

533 | 9 | pfleura2 | WRITE(*,*) " - Hybrid: geometries are described in zmat, but the gradient are used in cartesian" |

534 | 9 | pfleura2 | WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" |

535 | 9 | pfleura2 | WRITE(*,*) "" |

536 | 9 | pfleura2 | WRITE(*,*) " - RFO: Rational function optimization" |

537 | 9 | pfleura2 | WRITE(*,*) " - GDIIS: Geometry optimization using direct inversion in the iterative supspace" |

538 | 9 | pfleura2 | WRITE(*,*) "HUpdate: method to update the hessian. By default, it is Murtagh-Sargent Except for geom-" |

539 | 9 | pfleura2 | WRITE(*,*) " etry optimization where it is BFGS." |

540 | 9 | pfleura2 | WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization" |

541 | 9 | pfleura2 | WRITE(*,*) "" |

542 | 9 | pfleura2 | WRITE(*,*) "Smax: Maximum length of a step during path optimization" |

543 | 9 | pfleura2 | WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary" |

544 | 9 | pfleura2 | WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is" |

545 | 9 | pfleura2 | WRITE(*,*) " taken" |

546 | 9 | pfleura2 | WRITE(*,*) "" |

547 | 9 | pfleura2 | WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path" |

548 | 9 | pfleura2 | WRITE(*,*) " that is kept. FTan=1. corresponds to the full displacement, whereas FTan=0. gives a" |

549 | 9 | pfleura2 | WRITE(*,*) " displacement orthogonal to the path." |

550 | 9 | pfleura2 | WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of" |

551 | 9 | pfleura2 | WRITE(*,*) " reparameterization." |

552 | 9 | pfleura2 | WRITE(*,*) "IReparamT: When the path is not reparameterised at each iteration, this gives the frequency" |

553 | 9 | pfleura2 | WRITE(*,*) " of reparameterization of the tangents." |

554 | 9 | pfleura2 | WRITE(*,*) "" |

555 | 9 | pfleura2 | WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path. This option indicates" |

556 | 9 | pfleura2 | WRITE(*,*) " the first step where spline interpolation is used." |

557 | 9 | pfleura2 | WRITE(*,*) "BoxTol: Real between 0. and 1. When doing periodic calculations, it might happen that an" |

558 | 9 | pfleura2 | WRITE(*,*) " atom moves out of the unit cell. Path detects this by comparing the displacement to" |

559 | 9 | pfleura2 | WRITE(*,*) " boxtol: if an atom moves by more than Boxtol, then it is moved back into the unit cell." |

560 | 9 | pfleura2 | WRITE(*,*) " Default value: 0.5." |

561 | 9 | pfleura2 | WRITE(*,*) "FrozTol: (Real) This indicates the threshold to determine wether an atom moves between two" |

562 | 9 | pfleura2 | WRITE(*,*) " images. Default is 1e-4." |

563 | 9 | pfleura2 | WRITE(*,*) "" |

564 | 9 | pfleura2 | WRITE(*,*) "OptGeom: This INTEGER indicates the index of the geometry you want to optimize. If" |

565 | 9 | pfleura2 | WRITE(*,*) " OptGeom is set, then Opt'n Path performs a geometry optimization instead of a path" |

566 | 9 | pfleura2 | WRITE(*,*) " interpolation." |

567 | 9 | pfleura2 | WRITE(*,*) "GeomFile: Name of the file to print the geometries and their energy during a geometry opti-" |

568 | 9 | pfleura2 | WRITE(*,*) " mization. If this variable is not given then nothing is printed." |

569 | 9 | pfleura2 | WRITE(*,*) "AnaFile: Name of the file to print the values of the geometrical parameters that are monitored" |

570 | 9 | pfleura2 | WRITE(*,*) " if AnaGeom=.TRUE.. Default is PathName.dat" |

571 | 9 | pfleura2 | WRITE(*,*) "" |

572 | 9 | pfleura2 | WRITE(*,*) "GplotFile: Name of the gnuplot file to plot the evolution of the geometrical parameters that" |

573 | 9 | pfleura2 | WRITE(*,*) " are monitored if AnaGeom=.TRUE.. Default is PathName.gplot" |

574 | 9 | pfleura2 | WRITE(*,*) "" |

575 | 9 | pfleura2 | WRITE(*,*) "*******" |

576 | 9 | pfleura2 | WRITE(*,*) "* Arrays:" |

577 | 9 | pfleura2 | WRITE(*,*) "*******" |

578 | 9 | pfleura2 | WRITE(*,*) "" |

579 | 9 | pfleura2 | WRITE(*,*) "Rcov: Array containing the covalent radii of the first 86 elements. You can modify it using," |

580 | 9 | pfleura2 | WRITE(*,*) " rcov(6)=0.8." |

581 | 9 | pfleura2 | WRITE(*,*) "Mass: Array containing the atomic mass of the first 86 elements." |

582 | 9 | pfleura2 | WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations. If not given, Opt'n Path will" |

583 | 9 | pfleura2 | WRITE(*,*) " read the POTCAR file." |

584 | 9 | pfleura2 | WRITE(*,*) "" |

585 | 9 | pfleura2 | WRITE(*,*) "*******" |

586 | 9 | pfleura2 | WRITE(*,*) "* Flags:" |

587 | 9 | pfleura2 | WRITE(*,*) "*******" |

588 | 9 | pfleura2 | WRITE(*,*) "" |

589 | 9 | pfleura2 | WRITE(*,*) "MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE." |

590 | 9 | pfleura2 | WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE." |

591 | 9 | pfleura2 | WRITE(*,*) " unless for Coord equals CART." |

592 | 9 | pfleura2 | WRITE(*,*) "" |

593 | 9 | pfleura2 | WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products." |

594 | 9 | pfleura2 | WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants." |

595 | 9 | pfleura2 | WRITE(*,*) "CalcEProd: if TRUE the product energy will be computed. Default is FALSE. Not Compatible" |

596 | 9 | pfleura2 | WRITE(*,*) " with RunMode=Para." |

597 | 9 | pfleura2 | WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compat-" |

598 | 9 | pfleura2 | WRITE(*,*) " ible with RunMode=Para." |

599 | 9 | pfleura2 | WRITE(*,*) "" |

600 | 9 | pfleura2 | WRITE(*,*) "PathOnly: TRUE if one wants to generate the initial path, and stops." |

601 | 9 | pfleura2 | WRITE(*,*) "Align: If .FALSE., successive geometries along the path are not aligned on each other before" |

602 | 9 | pfleura2 | WRITE(*,*) " path interpolation. Default is .TRUE. if there are 4 atoms or more." |

603 | 9 | pfleura2 | WRITE(*,*) "Hinv: if True, then Hessian inversed is used." |

604 | 9 | pfleura2 | WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations." |

605 | 9 | pfleura2 | WRITE(*,*) "" |

606 | 9 | pfleura2 | WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of" |

607 | 9 | pfleura2 | WRITE(*,*) " the path." |

608 | 9 | pfleura2 | WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms. If True, a &frozenlist" |

609 | 9 | pfleura2 | WRITE(*,*) " namelist containing the list of frozen atoms must be given. If VASP is used, and frozen is" |

610 | 9 | pfleura2 | WRITE(*,*) " not given here, the program will use the F flags of the input geometry" |

611 | 9 | pfleura2 | WRITE(*,*) "FCart: True if one wants to describe some atoms using cartesian coordinates. *** Only used in" |

612 | 9 | pfleura2 | WRITE(*,*) " 'mixed' calculations. *** If True, a &cartlist namelist containing the list of cart atoms" |

613 | 9 | pfleura2 | WRITE(*,*) " must be given. By default, only frozen atoms are described in cartesian coordinates." |

614 | 9 | pfleura2 | WRITE(*,*) "" |

615 | 9 | pfleura2 | WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms." |

616 | 9 | pfleura2 | WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now." |

617 | 9 | pfleura2 | WRITE(*,*) "" |

618 | 9 | pfleura2 | WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files. See also" |

619 | 9 | pfleura2 | WRITE(*,*) " the POSCAR option. This can be used only if prog or input=VASP." |

620 | 9 | pfleura2 | WRITE(*,*) "AnaGeom: If TRUE, Opt'n Path will create a file .dat for geometries analysis. If True, " |

621 | 9 | pfleura2 | WRITE(*,*) "Opt'n Path will look for the AnaList namelist after the Path Namelist." |

622 | 9 | pfleura2 | WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry." |

623 | 9 | pfleura2 | WRITE(*,*) " If energy goes up, Smax = Smax*0.8, if not Smax = Smax*1.2. It is ensured that the" |

624 | 9 | pfleura2 | WRITE(*,*) " dynamical Smax is within [0.5*SMax_0 ,2*Smax_0 ]." |

625 | 9 | pfleura2 | WRITE(*,*) "" |

626 | 9 | pfleura2 | WRITE(*,*) "*******" |

627 | 9 | pfleura2 | WRITE(*,*) "* Additional namelists:" |

628 | 9 | pfleura2 | WRITE(*,*) "*******" |

629 | 9 | pfleura2 | WRITE(*,*) "" |

630 | 9 | pfleura2 | WRITE(*,*) "&cartlist list= ... &end This gives the list of atoms that are described using cartesian coor-" |

631 | 9 | pfleura2 | WRITE(*,*) " dinates. Read only if FCart=.TRUE.. To indicate an atom range, from 1 to 5 for example," |

632 | 9 | pfleura2 | WRITE(*,*) " one can put 1 -5 in the list. For example:" |

633 | 9 | pfleura2 | WRITE(*,*) " &cartlist list = 1 2 6 12 -20 &end" |

634 | 9 | pfleura2 | WRITE(*,*) " will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian." |

635 | 9 | pfleura2 | WRITE(*,*) "&Frozenlist list= ... &end This gives the list of atoms that are frozen during optimization." |

636 | 9 | pfleura2 | WRITE(*,*) " Read only if FFrozen=.TRUE.. To indicate an atom range, from 1 to 5 for example, one" |

637 | 9 | pfleura2 | WRITE(*,*) " can put 1 -5 in the list." |

638 | 9 | pfleura2 | WRITE(*,*) "" |

639 | 9 | pfleura2 | WRITE(*,*) "&Analist nb= ... &end list of variables for geometry analysis. If present and if Ana-" |

640 | 9 | pfleura2 | WRITE(*,*) " Geom=T then Opt'n Path will read it and print the values of the variable in a .dat file If" |

641 | 9 | pfleura2 | WRITE(*,*) " AnaGeom is T but Analist is not given, then Path.dat just contains the energy. Analist" |

642 | 9 | pfleura2 | WRITE(*,*) " contains the number of variables to monitor: Nb, and is followed by the description of the" |

643 | 9 | pfleura2 | WRITE(*,*) " variables among: b(ond) At1 At2 a(ngle) At1 At2 At3 d(ihedral) At1 At2 At3 At4 c NbAt" |

644 | 9 | pfleura2 | WRITE(*,*) " At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added at the" |

645 | 9 | pfleura2 | WRITE(*,*) " end of the real atoms of the system." |

646 | 9 | pfleura2 | !!!!!!!!!! end of included file |

647 | 1 | pfleura2 | WRITE(*,*) "" |

648 | 1 | pfleura2 | STOP |

649 | 1 | pfleura2 | END IF |

650 | 5 | pfleura2 | |

651 | 1 | pfleura2 | OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN') |

652 | 1 | pfleura2 | IOOUT=6 |

653 | 1 | pfleura2 | CASE (0) |

654 | 1 | pfleura2 | IOIN=5 |

655 | 1 | pfleura2 | IOOUT=6 |

656 | 1 | pfleura2 | END SELECT |

657 | 1 | pfleura2 | |

658 | 1 | pfleura2 | |

659 | 1 | pfleura2 | ! We initiliaze variables |

660 | 1 | pfleura2 | Pi=dacos(-1.0d0) |

661 | 1 | pfleura2 | Nat=1 |

662 | 1 | pfleura2 | Fact=1.1 |

663 | 1 | pfleura2 | IGeomRef=-1 |

664 | 1 | pfleura2 | NGeomI=1 |

665 | 1 | pfleura2 | NGeomF=8 |

666 | 1 | pfleura2 | IReparam=5 |

667 | 1 | pfleura2 | IReparamT=-1 |

668 | 1 | pfleura2 | ISpline=5 |

669 | 1 | pfleura2 | debug=.False. |

670 | 1 | pfleura2 | MW=.TRUE. |

671 | 1 | pfleura2 | bohr=.False. |

672 | 1 | pfleura2 | Coord='ZMAT' |

673 | 1 | pfleura2 | Step_method='RFO' |

674 | 1 | pfleura2 | DebugFile='Path.valid' |

675 | 1 | pfleura2 | PathName="Path" |

676 | 1 | pfleura2 | MaxCyc=50 |

677 | 1 | pfleura2 | SMax=0.1D0 |

678 | 1 | pfleura2 | SThresh=-1. |

679 | 1 | pfleura2 | GThresh=-1. |

680 | 1 | pfleura2 | FTan=0. |

681 | 1 | pfleura2 | Hinv=.TRUE. |

682 | 1 | pfleura2 | IniHUp=.FALSE. |

683 | 1 | pfleura2 | HupNeighbour=.FALSE. |

684 | 1 | pfleura2 | HesUpd="EMPTY" |

685 | 5 | pfleura2 | HUpdate="EMPTY" |

686 | 1 | pfleura2 | FFrozen=.FALSE. |

687 | 1 | pfleura2 | FCart=.FALSE. |

688 | 1 | pfleura2 | List=0 |

689 | 1 | pfleura2 | renum=.TRUE. |

690 | 1 | pfleura2 | OptReac=.FALSE. |

691 | 1 | pfleura2 | OptProd=.FALSE. |

692 | 1 | pfleura2 | CalcEReac=.FALSE. |

693 | 1 | pfleura2 | CalcEProd=.FALSE. |

694 | 1 | pfleura2 | EReac=0.d0 |

695 | 7 | pfleura2 | EProd=0.d0 |

696 | 1 | pfleura2 | OptGeom=-1 |

697 | 5 | pfleura2 | GeomFile="EMPTY" |

698 | 7 | pfleura2 | AnaGeom=.FALSE. |

699 | 7 | pfleura2 | AnaFile="EMPTY" |

700 | 8 | pfleura2 | GplotFile="EMPTY" |

701 | 7 | pfleura2 | Nb=0 |

702 | 7 | pfleura2 | NbCom=0 |

703 | 1 | pfleura2 | PathOnly=.FALSE. |

704 | 1 | pfleura2 | Prog='EMPTY' |

705 | 1 | pfleura2 | ProgExe='EMPTY' |

706 | 1 | pfleura2 | Runmode='SERIAL' |

707 | 1 | pfleura2 | Input='EMPTY' |

708 | 1 | pfleura2 | Poscar="POSCAR" |

709 | 1 | pfleura2 | AutoCart=.FALSE. |

710 | 1 | pfleura2 | VMD=.TRUE. |

711 | 1 | pfleura2 | WriteVASP=.FALSE. |

712 | 1 | pfleura2 | DynMaxStep=.FALSE. |

713 | 1 | pfleura2 | CalcName="EMPTY" |

714 | 1 | pfleura2 | ISuffix="EMPTY" |

715 | 1 | pfleura2 | OSuffix="EMPTY" |

716 | 1 | pfleura2 | NGintMax=10 |

717 | 1 | pfleura2 | Align=.TRUE. |

718 | 1 | pfleura2 | |

719 | 1 | pfleura2 | ! Number of point used to compute the length of the path, |

720 | 1 | pfleura2 | ! and to reparameterize the path |

721 | 1 | pfleura2 | NMaxPtPath=-1 |

722 | 1 | pfleura2 | FrozTol=1e-4 |

723 | 1 | pfleura2 | |

724 | 1 | pfleura2 | ! Read the namelist |

725 | 1 | pfleura2 | read (IOIN,path) |

726 | 1 | pfleura2 | |

727 | 5 | pfleura2 | debug=valid("Main").or.valid("Path") |

728 | 5 | pfleura2 | |

729 | 1 | pfleura2 | FCalcHess=valid("CalcHess") |

730 | 1 | pfleura2 | PathName=AdjustL(PathName) |

731 | 1 | pfleura2 | Coord=AdjustL(Coord) |

732 | 1 | pfleura2 | CALL UpCase(coord) |

733 | 1 | pfleura2 | |

734 | 1 | pfleura2 | Step_method=AdjustL(Step_method) |

735 | 1 | pfleura2 | CALL UpCase(Step_method) |

736 | 1 | pfleura2 | |

737 | 1 | pfleura2 | Prog=AdjustL(Prog) |

738 | 1 | pfleura2 | CALL UpCase(prog) |

739 | 1 | pfleura2 | |

740 | 1 | pfleura2 | Input=AdjustL(Input) |

741 | 1 | pfleura2 | CALL UpCase(input) |

742 | 1 | pfleura2 | |

743 | 1 | pfleura2 | Runmode=AdjustL(Runmode) |

744 | 1 | pfleura2 | Call UpCase(Runmode) |

745 | 1 | pfleura2 | IF (Runmode(1:4).EQ.'PARA') Runmode="PARA" |

746 | 1 | pfleura2 | |

747 | 1 | pfleura2 | if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN |

748 | 1 | pfleura2 | WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'" |

749 | 1 | pfleura2 | WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE" |

750 | 1 | pfleura2 | CalcEReac=.FALSE. |

751 | 1 | pfleura2 | CalcEProd=.FALSE. |

752 | 1 | pfleura2 | END IF |

753 | 1 | pfleura2 | |

754 | 1 | pfleura2 | If (IReparamT<=0) IReparamT=IReparam |

755 | 1 | pfleura2 | |

756 | 1 | pfleura2 | ! We take care of HesUpd flag |

757 | 1 | pfleura2 | ! this is needed because some users might still use it |

758 | 1 | pfleura2 | if (HesUpd.NE."EMPTY") THEN |

759 | 1 | pfleura2 | ! We are pityless: |

760 | 1 | pfleura2 | WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead" |

761 | 1 | pfleura2 | WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead" |

762 | 1 | pfleura2 | STOP |

763 | 1 | pfleura2 | END IF |

764 | 1 | pfleura2 | |

765 | 5 | pfleura2 | IF (HUpdate=='EMPTY') THEN |

766 | 5 | pfleura2 | IF (OptGeom>=1) THEN |

767 | 5 | pfleura2 | HupDate="BFGS" |

768 | 5 | pfleura2 | ELSE |

769 | 5 | pfleura2 | HUpdate="MS" |

770 | 5 | pfleura2 | END IF |

771 | 5 | pfleura2 | END IF |

772 | 1 | pfleura2 | |

773 | 7 | pfleura2 | If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN |

774 | 5 | pfleura2 | COORD='CART' |

775 | 5 | pfleura2 | END IF |

776 | 5 | pfleura2 | |

777 | 1 | pfleura2 | IF (COORD.EQ.'CART') THEN |

778 | 1 | pfleura2 | Renum=.FALSE. |

779 | 5 | pfleura2 | IF (OptGeom>0) THEN |

780 | 5 | pfleura2 | IGeomRef=OptGeom |

781 | 5 | pfleura2 | ELSE |

782 | 5 | pfleura2 | IGeomRef=1 |

783 | 5 | pfleura2 | END IF |

784 | 1 | pfleura2 | END IF |

785 | 1 | pfleura2 | |

786 | 5 | pfleura2 | |

787 | 1 | pfleura2 | if (Prog.EQ.'TEST2') Prog="CHAMFRE" |

788 | 5 | pfleura2 | if (Prog.EQ.'2D') Prog="TEST2D" |

789 | 3 | pfleura2 | if (Prog.EQ.'TEST3') Prog="LEPS" |

790 | 1 | pfleura2 | |

791 | 1 | pfleura2 | |

792 | 1 | pfleura2 | Select Case (Prog) |

793 | 5 | pfleura2 | CASE ("EMPTY","G09","GAUSSIAN") |

794 | 5 | pfleura2 | Prog='GAUSSIAN' |

795 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="g09" |

796 | 5 | pfleura2 | if (CalcName=="EMPTY") CalcName="Path" |

797 | 5 | pfleura2 | if (ISuffix=="EMPTY") ISuffix="com" |

798 | 5 | pfleura2 | if (OSuffix=="EMPTY") OSuffix="log" |

799 | 8 | pfleura2 | UnitProg="au" |

800 | 8 | pfleura2 | ConvE=au2kcal |

801 | 5 | pfleura2 | CASE ("G03") |

802 | 5 | pfleura2 | Prog='GAUSSIAN' |

803 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="g03" |

804 | 5 | pfleura2 | if (CalcName=="EMPTY") CalcName="Path" |

805 | 5 | pfleura2 | if (ISuffix=="EMPTY") ISuffix="com" |

806 | 5 | pfleura2 | if (OSuffix=="EMPTY") OSuffix="log" |

807 | 8 | pfleura2 | UnitProg="au" |

808 | 8 | pfleura2 | ConvE=au2kcal |

809 | 1 | pfleura2 | CASE ("EXT") |

810 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="./Ext.exe" |

811 | 1 | pfleura2 | if (CalcName=="EMPTY") CalcName="Ext" |

812 | 1 | pfleura2 | if (ISuffix=="EMPTY") ISuffix="in" |

813 | 1 | pfleura2 | if (OSuffix=="EMPTY") OSuffix="out" |

814 | 8 | pfleura2 | UnitProg="au" |

815 | 8 | pfleura2 | ConvE=au2kcal |

816 | 1 | pfleura2 | CASE ('MOPAC') |

817 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="mopac" |

818 | 1 | pfleura2 | if (CalcName=="EMPTY") CalcName="Path" |

819 | 1 | pfleura2 | if (ISuffix=="EMPTY") ISuffix="mop" |

820 | 1 | pfleura2 | if (OSuffix=="EMPTY") OSuffix="out" |

821 | 8 | pfleura2 | UnitProg="au" |

822 | 8 | pfleura2 | ConvE=au2kcal |

823 | 5 | pfleura2 | CASE ("SIESTA") |

824 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="siesta" |

825 | 1 | pfleura2 | if (CalcName=="EMPTY") CalcName="Path" |

826 | 5 | pfleura2 | if (ISuffix=="EMPTY") ISuffix="fdf" |

827 | 1 | pfleura2 | if (OSuffix=="EMPTY") OSuffix="out" |

828 | 8 | pfleura2 | UnitProg="eV" |

829 | 8 | pfleura2 | ConvE=eV2kcal |

830 | 5 | pfleura2 | CASE ("VASP") |

831 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="VASP" |

832 | 8 | pfleura2 | UnitProg="eV" |

833 | 8 | pfleura2 | ConvE=eV2kcal |

834 | 5 | pfleura2 | CASE ("TM","TURBOMOLE") |

835 | 5 | pfleura2 | Prog="TURBOMOLE" |

836 | 5 | pfleura2 | If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep" |

837 | 8 | pfleura2 | UnitProg="au" |

838 | 8 | pfleura2 | ConvE=au2kcal |

839 | 5 | pfleura2 | CASE ("TEST","CHAMFRE","TEST2D","LEPS") |

840 | 5 | pfleura2 | ProgExe="" |

841 | 8 | pfleura2 | UnitProg="au" |

842 | 8 | pfleura2 | ConvE=au2kcal |

843 | 5 | pfleura2 | CASE DEFAULT |

844 | 5 | pfleura2 | WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP" |

845 | 5 | pfleura2 | STOP |

846 | 1 | pfleura2 | END SELECT |

847 | 1 | pfleura2 | |

848 | 1 | pfleura2 | if (Input.EQ.'EMPTY') THEN |

849 | 1 | pfleura2 | Select Case (Prog) |

850 | 1 | pfleura2 | CASE ("VASP") |

851 | 1 | pfleura2 | Input=Prog |

852 | 1 | pfleura2 | CASE ("CHAMFRE") |

853 | 1 | pfleura2 | Input="VASP" |

854 | 1 | pfleura2 | CASE DEFAULT |

855 | 1 | pfleura2 | Input='XYZ' |

856 | 1 | pfleura2 | END SELECT |

857 | 1 | pfleura2 | WRITE(*,*) "Input has been set to the default: ",INPUT |

858 | 1 | pfleura2 | END IF |

859 | 1 | pfleura2 | |

860 | 1 | pfleura2 | WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) |

861 | 1 | pfleura2 | |

862 | 1 | pfleura2 | IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN |

863 | 1 | pfleura2 | WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode." |

864 | 1 | pfleura2 | STOP |

865 | 1 | pfleura2 | END IF |

866 | 1 | pfleura2 | |

867 | 1 | pfleura2 | if (OptGeom.GE.1) THEN |

868 | 5 | pfleura2 | FPrintGeom=.FALSE. |

869 | 5 | pfleura2 | If (GeomFile/='EMPTY') THEN |

870 | 5 | pfleura2 | FPrintGeom=.TRUE. |

871 | 5 | pfleura2 | END IF |

872 | 1 | pfleura2 | IF (IGeomRef.LE.0) THEN |

873 | 1 | pfleura2 | IGeomRef=OptGeom |

874 | 1 | pfleura2 | ELSE |

875 | 1 | pfleura2 | IF (IGeomRef.NE.OptGeom) THEN |

876 | 1 | pfleura2 | WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !" |

877 | 1 | pfleura2 | WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom |

878 | 1 | pfleura2 | |

879 | 1 | pfleura2 | WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !" |

880 | 1 | pfleura2 | WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom |

881 | 1 | pfleura2 | STOP |

882 | 1 | pfleura2 | END IF |

883 | 1 | pfleura2 | END IF |

884 | 1 | pfleura2 | END IF |

885 | 1 | pfleura2 | |

886 | 5 | pfleura2 | ALLOCATE(Cart(Nat)) |

887 | 5 | pfleura2 | Cart=0 |

888 | 1 | pfleura2 | |

889 | 5 | pfleura2 | ALLOCATE(Frozen(Nat)) |

890 | 5 | pfleura2 | Frozen=0 |

891 | 1 | pfleura2 | |

892 | 5 | pfleura2 | IF (FCart) THEN |

893 | 5 | pfleura2 | ! We rewind the file to be sure to browse all of it to read the Cart namelist |

894 | 5 | pfleura2 | REWIND(IOIN) |

895 | 5 | pfleura2 | List=0 |

896 | 5 | pfleura2 | READ(IOIN,cartlist) |

897 | 5 | pfleura2 | IF (COORD.NE.'CART') THEN |

898 | 5 | pfleura2 | Cart=List(1:Nat) |

899 | 5 | pfleura2 | if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart |

900 | 1 | pfleura2 | ELSE |

901 | 5 | pfleura2 | WRITE(*,*) "FCart=T is not allowed when Coord='CART'" |

902 | 5 | pfleura2 | WRITE(*,*) "I will discard the cartlist" |

903 | 1 | pfleura2 | Cart=0 |

904 | 1 | pfleura2 | END IF |

905 | 1 | pfleura2 | END IF |

906 | 1 | pfleura2 | |

907 | 5 | pfleura2 | if (FFrozen) THEN |

908 | 1 | pfleura2 | |

909 | 5 | pfleura2 | REWIND(IOIN) |

910 | 5 | pfleura2 | List=0 |

911 | 5 | pfleura2 | READ(IOIN,Frozenlist) |

912 | 5 | pfleura2 | Frozen=List(1:Nat) |

913 | 5 | pfleura2 | if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen |

914 | 5 | pfleura2 | END IF |

915 | 5 | pfleura2 | |

916 | 1 | pfleura2 | If (FCart) Call Expand(Cart,NCart,Nat) |

917 | 1 | pfleura2 | IF (FFrozen) Call Expand(Frozen,NFroz,Nat) |

918 | 1 | pfleura2 | |

919 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart |

920 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen |

921 | 1 | pfleura2 | |

922 | 7 | pfleura2 | if (AnaGeom) THEN |

923 | 7 | pfleura2 | REWIND(IOIN) |

924 | 7 | pfleura2 | READ(IOIN,AnaList,IOSTAT=IoS) |

925 | 7 | pfleura2 | IF (IOS==0) THEN ! we have a AnaList namelist |

926 | 7 | pfleura2 | Call ReadAnaList |

927 | 8 | pfleura2 | ELSE |

928 | 8 | pfleura2 | ! No AnaList namelist, we will print only the energy |

929 | 8 | pfleura2 | FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)' |

930 | 7 | pfleura2 | END IF |

931 | 7 | pfleura2 | IF (AnaFile=="EMPTY") THEN |

932 | 7 | pfleura2 | AnaFile=TRIM(PathName) // '.dat' |

933 | 7 | pfleura2 | END IF |

934 | 7 | pfleura2 | OPEN(IODat,File=AnaFile) |

935 | 8 | pfleura2 | IF (GplotFile=="EMPTY") THEN |

936 | 8 | pfleura2 | GplotFile=TRIM(PathName) // '.gplot' |

937 | 8 | pfleura2 | END IF |

938 | 8 | pfleura2 | |

939 | 8 | pfleura2 | OPEN(IODat,File=AnaFile) |

940 | 8 | pfleura2 | If (OptGeom<=0) THEN |

941 | 8 | pfleura2 | OPEN(IOGplot,File=GPlotFile) |

942 | 8 | pfleura2 | WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist" |

943 | 8 | pfleura2 | WRITE(IoGplot,'(A)') " set pointsize 2" |

944 | 8 | pfleura2 | WRITE(IoGplot,'(A)') " xcol=1" |

945 | 8 | pfleura2 | WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3 |

946 | 8 | pfleura2 | END IF |

947 | 7 | pfleura2 | END IF |

948 | 7 | pfleura2 | |

949 | 7 | pfleura2 | |

950 | 7 | pfleura2 | |

951 | 1 | pfleura2 | if (NMaxPtPath.EQ.-1) then |

952 | 1 | pfleura2 | NMaxPtPath=min(50*NGeomF,2000) |

953 | 1 | pfleura2 | NMaxPtPath=max(20*NGeomF,NMaxPtPath) |

954 | 1 | pfleura2 | end if |

955 | 1 | pfleura2 | |

956 | 1 | pfleura2 | ! NMaxPtPath=10000 |

957 | 1 | pfleura2 | |

958 | 1 | pfleura2 | ! We ensure that FTan is in [0.:1.] |

959 | 1 | pfleura2 | FTan=min(abs(FTan),1.d0) |

960 | 1 | pfleura2 | |

961 | 1 | pfleura2 | ! PFL 2011 Mar 14 -> |

962 | 1 | pfleura2 | ! Added some printing for analyses with Anapath |

963 | 1 | pfleura2 | if (debug) THEN |

964 | 1 | pfleura2 | WRITE(IOOUT,path) |

965 | 1 | pfleura2 | ELSE |

966 | 1 | pfleura2 | ! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC |

967 | 1 | pfleura2 | ! and PATHNAME |

968 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat |

969 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI |

970 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF |

971 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC |

972 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME) |

973 | 1 | pfleura2 | END IF |

974 | 1 | pfleura2 | |

975 | 1 | pfleura2 | FTan=1.-FTan |

976 | 1 | pfleura2 | |

977 | 1 | pfleura2 | !Thresholds... |

978 | 1 | pfleura2 | if (SThresh.LE.0) SThresh=0.01d0 |

979 | 1 | pfleura2 | If (GThresh.LE.0) GThresh=SThresh/4. |

980 | 1 | pfleura2 | SThreshRms=Sthresh*2./3. |

981 | 1 | pfleura2 | GThreshRms=Gthresh*2./3. |

982 | 1 | pfleura2 | |

983 | 1 | pfleura2 | ! First batch of allocations. Arrays that depend only on NGeomI, Nat |

984 | 1 | pfleura2 | ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) |

985 | 1 | pfleura2 | ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat)) |

986 | 1 | pfleura2 | ALLOCATE(MassAt(NAt),OrderInv(Nat)) |

987 | 1 | pfleura2 | |

988 | 7 | pfleura2 | AtName="" |

989 | 7 | pfleura2 | MassAt=1. |

990 | 7 | pfleura2 | |

991 | 8 | pfleura2 | ! As we have used rewind many times, the position in the input file |

992 | 8 | pfleura2 | ! is blurry (at least !) we thus have to go back to the Geometries description |

993 | 8 | pfleura2 | ! (TO DO: use a better Input parser !) |

994 | 8 | pfleura2 | |

995 | 8 | pfleura2 | REWIND(IOIN) |

996 | 8 | pfleura2 | ! We search the '&path' namelist |

997 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

998 | 8 | pfleura2 | Fini=.FALSE. |

999 | 8 | pfleura2 | DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0) |

1000 | 8 | pfleura2 | if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line) |

1001 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

1002 | 8 | pfleura2 | END DO |

1003 | 8 | pfleura2 | if (debug) WRITE(*,*) "Path found:",TRIM(LINE) |

1004 | 8 | pfleura2 | ! We are on the &path line, we move to the end |

1005 | 8 | pfleura2 | Call NoComment(Line) |

1006 | 8 | pfleura2 | DO WHILE ((Ios==0).AND.InString(Line,"/")==0) |

1007 | 8 | pfleura2 | if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line) |

1008 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

1009 | 8 | pfleura2 | Call NoComment(Line) |

1010 | 8 | pfleura2 | Call noString(Line) |

1011 | 8 | pfleura2 | END DO |

1012 | 8 | pfleura2 | ! We are at then end of &Path / namelist |

1013 | 8 | pfleura2 | ! We scan for other Namelists |

1014 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

1015 | 8 | pfleura2 | if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line) |

1016 | 8 | pfleura2 | DO WHILE ((IoS==0).AND.(Index(Line,'&')/=0)) |

1017 | 8 | pfleura2 | Call NoComment(Line) |

1018 | 8 | pfleura2 | if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line) |

1019 | 8 | pfleura2 | Idx=InString(Line,'Analist') |

1020 | 8 | pfleura2 | DO WHILE ((Ios==0).AND.InString(Line,"/")==0) |

1021 | 8 | pfleura2 | if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line) |

1022 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

1023 | 8 | pfleura2 | Call NoComment(Line) |

1024 | 8 | pfleura2 | END DO |

1025 | 8 | pfleura2 | If (Idx/=0) THEN |

1026 | 8 | pfleura2 | ! we have just read &Analist namelist, we have to skip the variables description |

1027 | 8 | pfleura2 | DO I=1,Nb |

1028 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

1029 | 8 | pfleura2 | if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE) |

1030 | 8 | pfleura2 | END DO |

1031 | 8 | pfleura2 | END IF |

1032 | 8 | pfleura2 | READ(IOIN,'(A)',Iostat=Ios) Line |

1033 | 8 | pfleura2 | if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE) |

1034 | 8 | pfleura2 | END DO |

1035 | 8 | pfleura2 | BACKSPACE(IOIN) |

1036 | 8 | pfleura2 | |

1037 | 1 | pfleura2 | ! We read the initial geometries |

1038 | 1 | pfleura2 | Call Read_Geom(input) |

1039 | 1 | pfleura2 | |

1040 | 1 | pfleura2 | ! We convert atome names into atom mass number |

1041 | 1 | pfleura2 | DO I=1,NAt |

1042 | 1 | pfleura2 | ! WRITE(*,*) 'I,AtName',I,AtName(I) |

1043 | 1 | pfleura2 | Atome(I)=ConvertNumAt(AtName(I)) |

1044 | 1 | pfleura2 | END DO |

1045 | 1 | pfleura2 | |

1046 | 7 | pfleura2 | If (AnaGeom) THEN |

1047 | 8 | pfleura2 | If (NbVar>0) THEN |

1048 | 7 | pfleura2 | ! We print the description of the varialbes in AnaFile |

1049 | 8 | pfleura2 | Call PrintAnaList(IoDat) |

1050 | 8 | pfleura2 | If (OptGeom<=0) THEN |

1051 | 8 | pfleura2 | WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")" |

1052 | 8 | pfleura2 | ELSE |

1053 | 8 | pfleura2 | WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")" |

1054 | 8 | pfleura2 | END IF |

1055 | 8 | pfleura2 | ELSE |

1056 | 8 | pfleura2 | If (OptGeom<=0) THEN |

1057 | 8 | pfleura2 | WRITE(IoDat,'(A)') "# s Erel(kcal/mol) E (" // TRIM(UnitProg) //")" |

1058 | 8 | pfleura2 | ELSE |

1059 | 8 | pfleura2 | WRITE(IoDat,'(A)') "# Iter Erel(kcal/mol) E (" // TRIM(UnitProg) //")" |

1060 | 8 | pfleura2 | END IF |

1061 | 8 | pfleura2 | END IF |

1062 | 7 | pfleura2 | END IF |

1063 | 1 | pfleura2 | |

1064 | 1 | pfleura2 | ! PFL 23 Sep 2008 -> |

1065 | 1 | pfleura2 | ! We check that frozen atoms are really frozen, ie that their coordinates do not change |

1066 | 1 | pfleura2 | ! between the first geometry and the others. |

1067 | 1 | pfleura2 | IF (NFroz.GT.0) THEN |

1068 | 1 | pfleura2 | ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat)) |

1069 | 1 | pfleura2 | ALLOCATE(x2(nat),y2(nat),z2(nat)) |

1070 | 1 | pfleura2 | ListAtFroz=.FALSE. |

1071 | 1 | pfleura2 | |

1072 | 1 | pfleura2 | x1(1:Nat)=XyzgeomI(1,1,1:Nat) |

1073 | 1 | pfleura2 | y1(1:Nat)=XyzgeomI(1,2,1:Nat) |

1074 | 1 | pfleura2 | z1(1:Nat)=XyzgeomI(1,3,1:Nat) |

1075 | 1 | pfleura2 | |

1076 | 1 | pfleura2 | SELECT CASE (NFroz) |

1077 | 1 | pfleura2 | ! We have Frozen atoms |

1078 | 1 | pfleura2 | ! PFL 28 Oct 2008 -> |

1079 | 1 | pfleura2 | ! It might happen that the geometries are translated/rotated. |

1080 | 1 | pfleura2 | ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first |

1081 | 1 | pfleura2 | ! three frozen atoms. |

1082 | 1 | pfleura2 | |

1083 | 1 | pfleura2 | CASE (3:) |

1084 | 1 | pfleura2 | DO I=1,NFroz |

1085 | 1 | pfleura2 | ListAtFroz(Frozen(I))=.TRUE. |

1086 | 1 | pfleura2 | END DO |

1087 | 1 | pfleura2 | CASE (2) |

1088 | 1 | pfleura2 | IAt1=Frozen(1) |

1089 | 1 | pfleura2 | IAt2=Frozen(2) |

1090 | 1 | pfleura2 | ListAtFroz(Iat1)=.TRUE. |

1091 | 1 | pfleura2 | ListAtFroz(Iat2)=.TRUE. |

1092 | 1 | pfleura2 | x2(1)=x1(Iat2)-x1(Iat1) |

1093 | 1 | pfleura2 | y2(1)=y1(Iat2)-y1(Iat1) |

1094 | 1 | pfleura2 | z2(1)=z1(Iat2)-z1(Iat1) |

1095 | 1 | pfleura2 | Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) |

1096 | 1 | pfleura2 | Dist=-1. |

1097 | 1 | pfleura2 | ! We scan all atoms to find one that is not aligned with Iat1-Iat2 |

1098 | 1 | pfleura2 | DO I=1,Nat |

1099 | 1 | pfleura2 | if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle |

1100 | 1 | pfleura2 | x2(2)=x1(Iat2)-x1(I) |

1101 | 1 | pfleura2 | y2(2)=y1(Iat2)-y1(I) |

1102 | 1 | pfleura2 | z2(2)=z1(Iat2)-z1(I) |

1103 | 1 | pfleura2 | Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2) |

1104 | 1 | pfleura2 | ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2) |

1105 | 1 | pfleura2 | if (abs(ca)<0.9) THEN |

1106 | 1 | pfleura2 | IF (Norm2>Dist) THEN |

1107 | 1 | pfleura2 | Iat3=I |

1108 | 1 | pfleura2 | Dist=Norm2 |

1109 | 1 | pfleura2 | END IF |

1110 | 1 | pfleura2 | END IF |

1111 | 1 | pfleura2 | END DO |

1112 | 1 | pfleura2 | ListAtFroz(Iat3)=.TRUE. |

1113 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3 |

1114 | 1 | pfleura2 | CASE (1) |

1115 | 1 | pfleura2 | IAt1=Frozen(1) |

1116 | 1 | pfleura2 | ListAtFroz(Iat1)=.TRUE. |

1117 | 1 | pfleura2 | Dist=-1. |

1118 | 1 | pfleura2 | ! We scan all atoms to find one that is further from At1 |

1119 | 1 | pfleura2 | DO I=1,Nat |

1120 | 1 | pfleura2 | if (I.EQ.Iat1) Cycle |

1121 | 1 | pfleura2 | x2(1)=x1(Iat1)-x1(I) |

1122 | 1 | pfleura2 | y2(1)=y1(Iat1)-y1(I) |

1123 | 1 | pfleura2 | z2(1)=z1(Iat1)-z1(I) |

1124 | 1 | pfleura2 | Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) |

1125 | 1 | pfleura2 | IF (Norm1>Dist) THEN |

1126 | 1 | pfleura2 | Iat2=I |

1127 | 1 | pfleura2 | Dist=Norm1 |

1128 | 1 | pfleura2 | END IF |

1129 | 1 | pfleura2 | END DO |

1130 | 1 | pfleura2 | ListAtFroz(Iat2)=.TRUE. |

1131 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2 |

1132 | 1 | pfleura2 | x2(1)=x1(Iat2)-x1(Iat1) |

1133 | 1 | pfleura2 | y2(1)=y1(Iat2)-y1(Iat1) |

1134 | 1 | pfleura2 | z2(1)=z1(Iat2)-z1(Iat1) |

1135 | 1 | pfleura2 | Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2) |

1136 | 1 | pfleura2 | Dist=-1. |

1137 | 1 | pfleura2 | ! We scan all atoms to find one that is not aligned with Iat1-Iat2 |

1138 | 8 | pfleura2 | Iat3=1 |

1139 | 8 | pfleura2 | DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat)) |

1140 | 8 | pfleura2 | Iat3=Iat3+1 |

1141 | 8 | pfleura2 | END DO |

1142 | 1 | pfleura2 | DO I=1,Nat |

1143 | 1 | pfleura2 | if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle |

1144 | 1 | pfleura2 | x2(2)=x1(Iat2)-x1(I) |

1145 | 1 | pfleura2 | y2(2)=y1(Iat2)-y1(I) |

1146 | 1 | pfleura2 | z2(2)=z1(Iat2)-z1(I) |

1147 | 1 | pfleura2 | Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2) |

1148 | 1 | pfleura2 | ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2) |

1149 | 1 | pfleura2 | if (abs(ca)<0.9) THEN |

1150 | 1 | pfleura2 | IF (Norm2>Dist) THEN |

1151 | 1 | pfleura2 | Iat3=I |

1152 | 1 | pfleura2 | Dist=Norm2 |

1153 | 1 | pfleura2 | END IF |

1154 | 1 | pfleura2 | END IF |

1155 | 1 | pfleura2 | END DO |

1156 | 1 | pfleura2 | ListAtFroz(Iat3)=.TRUE. |

1157 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3 |

1158 | 1 | pfleura2 | END SELECT |

1159 | 1 | pfleura2 | |

1160 | 1 | pfleura2 | DO I=2,NGeomI |

1161 | 1 | pfleura2 | ! First, we align the Ith geometry on the first one, using only |

1162 | 1 | pfleura2 | ! the positions of the "frozen" atoms. (see above: it might be that |

1163 | 8 | pfleura2 | ! ListAtFroz contains non frozen atoms used to align the geometries |

1164 | 1 | pfleura2 | x2(1:Nat)=XyzgeomI(I,1,1:Nat) |

1165 | 1 | pfleura2 | y2(1:Nat)=XyzgeomI(I,2,1:Nat) |

1166 | 1 | pfleura2 | z2(1:Nat)=XyzgeomI(I,3,1:Nat) |

1167 | 1 | pfleura2 | Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot) |

1168 | 1 | pfleura2 | |

1169 | 1 | pfleura2 | FChkFrozen=.FALSE. |

1170 | 1 | pfleura2 | DO J=1,NFroz |

1171 | 1 | pfleura2 | IAt=Frozen(J) |

1172 | 1 | pfleura2 | IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. & |

1173 | 1 | pfleura2 | (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR. & |

1174 | 1 | pfleura2 | (abs(z1(Iat)-z2(Iat)).GT.FrozTol)) THEN |

1175 | 1 | pfleura2 | IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****" |

1176 | 1 | pfleura2 | FChkFrozen=.TRUE. |

1177 | 1 | pfleura2 | WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change" |

1178 | 1 | pfleura2 | WRITE(*,*) "from geometry 1 to geometry ",I |

1179 | 1 | pfleura2 | END IF |

1180 | 1 | pfleura2 | END DO |

1181 | 1 | pfleura2 | END DO |

1182 | 1 | pfleura2 | IF (FChkFrozen) THEN |

1183 | 1 | pfleura2 | WRITE(*,*) "Some frozen atoms are not frozen... check this and play again" |

1184 | 1 | pfleura2 | STOP |

1185 | 1 | pfleura2 | END IF |

1186 | 1 | pfleura2 | END IF |

1187 | 1 | pfleura2 | |

1188 | 1 | pfleura2 | |

1189 | 1 | pfleura2 | N3at=3*Nat |

1190 | 1 | pfleura2 | NFree=N3at-6 ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES. |

1191 | 1 | pfleura2 | |

1192 | 1 | pfleura2 | |

1193 | 5 | pfleura2 | Call ReadInput |

1194 | 1 | pfleura2 | |

1195 | 1 | pfleura2 | IF (COORD=="MIXED") Call TestCart(AutoCart) |

1196 | 1 | pfleura2 | |

1197 | 1 | pfleura2 | SELECT CASE(NFroz) |

1198 | 1 | pfleura2 | CASE (2) |

1199 | 1 | pfleura2 | IntFroz=1 |

1200 | 1 | pfleura2 | CASE (3) |

1201 | 1 | pfleura2 | IntFroz=3 |

1202 | 1 | pfleura2 | CASE (4:) |

1203 | 1 | pfleura2 | IntFroz=(NFroz-2)*3 !(NFroz-3)*3+3 |

1204 | 1 | pfleura2 | END SELECT |

1205 | 1 | pfleura2 | |

1206 | 1 | pfleura2 | if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "& |

1207 | 1 | pfleura2 | ,TRIM(COORD),Ncart,NCoord,N3at,NFroz |

1208 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart |

1209 | 1 | pfleura2 | |

1210 | 1 | pfleura2 | if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN |

1211 | 1 | pfleura2 | WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT." |

1212 | 1 | pfleura2 | COORD="ZMAT" |

1213 | 1 | pfleura2 | END IF |

1214 | 1 | pfleura2 | |

1215 | 1 | pfleura2 | IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz |

1216 | 1 | pfleura2 | |

1217 | 1 | pfleura2 | SELECT CASE(COORD) |

1218 | 1 | pfleura2 | CASE ('ZMAT') |

1219 | 1 | pfleura2 | NCoord=Nfree |

1220 | 1 | pfleura2 | ALLOCATE(dzdc(3,nat,3,nat)) |

1221 | 1 | pfleura2 | CASE ('MIXED') |

1222 | 1 | pfleura2 | SELECT CASE (NCart+NFroz) |

1223 | 1 | pfleura2 | CASE (1) |

1224 | 1 | pfleura2 | NCoord=N3At-3 |

1225 | 1 | pfleura2 | CASE (2) |

1226 | 1 | pfleura2 | NCoord=N3At-1 |

1227 | 1 | pfleura2 | CASE (3:) |

1228 | 1 | pfleura2 | NCoord=N3At |

1229 | 1 | pfleura2 | END SELECT |

1230 | 1 | pfleura2 | ALLOCATE(dzdc(3,nat,3,nat)) |

1231 | 1 | pfleura2 | CASE ('BAKER') |

1232 | 1 | pfleura2 | Symmetry_elimination=0 |

1233 | 1 | pfleura2 | NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6 |

1234 | 1 | pfleura2 | ALLOCATE(BMat_BakerT(3*Nat,NCoord)) |

1235 | 1 | pfleura2 | ALLOCATE(BTransInv(NCoord,3*Nat)) |

1236 | 1 | pfleura2 | ALLOCATE(BTransInv_local(NCoord,3*Nat)) |

1237 | 1 | pfleura2 | CASE ('HYBRID') |

1238 | 1 | pfleura2 | NCoord=N3at |

1239 | 1 | pfleura2 | CASE ('CART') |

1240 | 1 | pfleura2 | NCoord=N3at |

1241 | 1 | pfleura2 | END SELECT |

1242 | 1 | pfleura2 | |

1243 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Path, L826, Coord, NCart,NCoord, N3At",TRIM(COORD),Ncart, NCoord,N3at |

1244 | 1 | pfleura2 | |

1245 | 1 | pfleura2 | ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord)) |

1246 | 1 | pfleura2 | ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat)) |

1247 | 1 | pfleura2 | ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord)) |

1248 | 1 | pfleura2 | ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in Path_module.f90 |

1249 | 1 | pfleura2 | ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord)) |

1250 | 1 | pfleura2 | ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord)) |

1251 | 1 | pfleura2 | ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in Path_module.f90 |

1252 | 1 | pfleura2 | ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF)) |

1253 | 1 | pfleura2 | |

1254 | 1 | pfleura2 | ZeroVec=0.d0 |

1255 | 1 | pfleura2 | DO I =1, NGeomF |

1256 | 1 | pfleura2 | IntTangent(I,:)=0.d0 |

1257 | 1 | pfleura2 | END DO |

1258 | 1 | pfleura2 | |

1259 | 1 | pfleura2 | if (debug) THEN |

1260 | 1 | pfleura2 | Print *, 'L886, IntTangent(1,:)=' |

1261 | 1 | pfleura2 | Print *, IntTangent(1,:) |

1262 | 1 | pfleura2 | END IF |

1263 | 1 | pfleura2 | |

1264 | 1 | pfleura2 | if (.NOT.OptReac) Energies(1)=Ereac |

1265 | 1 | pfleura2 | if (.NOT.OptProd) Energies(NGeomF)=EProd |

1266 | 1 | pfleura2 | MaxStep=SMax |

1267 | 1 | pfleura2 | |

1268 | 1 | pfleura2 | ! We initialize the mass array |

1269 | 1 | pfleura2 | if (MW) THEN |

1270 | 1 | pfleura2 | WRITE(*,*) 'Working in MW coordinates' |

1271 | 1 | pfleura2 | DO I=1,Nat |

1272 | 1 | pfleura2 | massAt(I)=Mass(Atome(I)) |

1273 | 1 | pfleura2 | END DO |

1274 | 1 | pfleura2 | ELSE |

1275 | 1 | pfleura2 | DO I=1,Nat |

1276 | 1 | pfleura2 | MassAt(I)=1.d0 |

1277 | 1 | pfleura2 | END DO |

1278 | 1 | pfleura2 | END IF |

1279 | 1 | pfleura2 | |

1280 | 1 | pfleura2 | WRITE(*,*) "Prog=",TRIM(PROG) |

1281 | 1 | pfleura2 | |

1282 | 1 | pfleura2 | ALLOCATE(FrozAtoms(Nat)) |

1283 | 1 | pfleura2 | |

1284 | 1 | pfleura2 | ! IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz |

1285 | 1 | pfleura2 | ! IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen |

1286 | 1 | pfleura2 | ! IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms |

1287 | 1 | pfleura2 | IF (NFroz.EQ.0) THEN |

1288 | 1 | pfleura2 | FrozAtoms=.TRUE. |

1289 | 1 | pfleura2 | ELSE |

1290 | 1 | pfleura2 | I=1 |

1291 | 1 | pfleura2 | NFr=0 |

1292 | 1 | pfleura2 | FrozAtoms=.False. |

1293 | 1 | pfleura2 | DO WHILE (Frozen(I).NE.0) |

1294 | 1 | pfleura2 | IF (Frozen(I).LT.0) THEN |

1295 | 1 | pfleura2 | DO J=Frozen(I-1),abs(Frozen(I)) |

1296 | 1 | pfleura2 | IF (.NOT.FrozAtoms(J)) THEN |

1297 | 1 | pfleura2 | NFr=NFr+1 |

1298 | 1 | pfleura2 | FrozAtoms(J)=.TRUE. |

1299 | 1 | pfleura2 | END IF |

1300 | 1 | pfleura2 | END DO |

1301 | 1 | pfleura2 | ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN |

1302 | 1 | pfleura2 | FrozAtoms(Frozen(I))=.TRUE. |

1303 | 1 | pfleura2 | NFr=NFr+1 |

1304 | 1 | pfleura2 | END IF |

1305 | 1 | pfleura2 | I=I+1 |

1306 | 1 | pfleura2 | END DO |

1307 | 1 | pfleura2 | IF (NFr.NE.NFroz) THEN |

1308 | 1 | pfleura2 | WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr |

1309 | 1 | pfleura2 | STOP |

1310 | 1 | pfleura2 | END IF |

1311 | 1 | pfleura2 | END IF |

1312 | 1 | pfleura2 | |

1313 | 5 | pfleura2 | ! if (debug) THEN |

1314 | 1 | pfleura2 | !Print *, 'L968, IntTangent(1,:)=' |

1315 | 1 | pfleura2 | !Print *, IntTangent(1,:) |

1316 | 5 | pfleura2 | ! END IF |

1317 | 1 | pfleura2 | |

1318 | 1 | pfleura2 | ! We have read everything we needed... time to work :-) |

1319 | 1 | pfleura2 | IOpt=0 |

1320 | 1 | pfleura2 | FirstTimePathCreate = .TRUE. ! Initialized. |

1321 | 1 | pfleura2 | Call PathCreate(IOpt) ! IntTangent is also computed in PathCreate. |

1322 | 1 | pfleura2 | FirstTimePathCreate = .FALSE. |

1323 | 1 | pfleura2 | |

1324 | 1 | pfleura2 | if (debug) THEN |

1325 | 1 | pfleura2 | Print *, 'L980, IntTangent(1,:)=' |

1326 | 1 | pfleura2 | Print *, IntTangent(1,:) |

1327 | 1 | pfleura2 | END IF |

1328 | 1 | pfleura2 | |

1329 | 1 | pfleura2 | ! PFL 30 Mar 2008 -> |

1330 | 1 | pfleura2 | ! The following is now allocated in Calc_Baker. This is more logical to me |

1331 | 1 | pfleura2 | ! IF (COORD .EQ. "BAKER") Then |

1332 | 1 | pfleura2 | ! ALLOCATE(XprimitiveF(NgeomF,NbCoords)) |

1333 | 1 | pfleura2 | ! ALLOCATE(UMatF(NGeomF,NbCoords,NCoord)) |

1334 | 1 | pfleura2 | ! ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat)) |

1335 | 1 | pfleura2 | ! END IF |

1336 | 1 | pfleura2 | ! <- PFL 30 mar 2008 |

1337 | 1 | pfleura2 | |

1338 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1339 | 1 | pfleura2 | ! |

1340 | 1 | pfleura2 | ! |

1341 | 1 | pfleura2 | ! OptGeom is the index of the geometry which is to be optimized. |

1342 | 1 | pfleura2 | IF (Optgeom.GT.0) THEN |

1343 | 1 | pfleura2 | Flag_Opt_Geom=.TRUE. |

1344 | 1 | pfleura2 | SELECT CASE(Coord) |

1345 | 1 | pfleura2 | CASE ('CART','HYBRID') |

1346 | 1 | pfleura2 | !!! CAUTION : PFL 29.AUG.2008 -> |

1347 | 1 | pfleura2 | ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3) |

1348 | 1 | pfleura2 | ! This should be made consistent !!!!!!!!!!!!!!! |

1349 | 5 | pfleura2 | GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) |

1350 | 1 | pfleura2 | ! GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) |

1351 | 1 | pfleura2 | !!! <- CAUTION : PFL 29.AUG.2008 |

1352 | 1 | pfleura2 | CASE ('ZMAT','MIXED') |

1353 | 1 | pfleura2 | !GeomTmp=IntCoordF(OptGeom,:) |

1354 | 1 | pfleura2 | GeomTmp=IntCoordI(OptGeom,:) |

1355 | 1 | pfleura2 | CASE ('BAKER') |

1356 | 1 | pfleura2 | !GeomTmp=IntCoordF(OptGeom,:) |

1357 | 1 | pfleura2 | GeomTmp=IntCoordI(OptGeom,:) |

1358 | 1 | pfleura2 | CASE DEFAULT |

1359 | 1 | pfleura2 | WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP" |

1360 | 1 | pfleura2 | STOP |

1361 | 1 | pfleura2 | END SELECT |

1362 | 1 | pfleura2 | |

1363 | 1 | pfleura2 | ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker |

1364 | 1 | pfleura2 | Flag_Opt_Geom = .TRUE. |

1365 | 1 | pfleura2 | Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method) |

1366 | 1 | pfleura2 | |

1367 | 1 | pfleura2 | WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom |

1368 | 9 | pfleura2 | if (debug) WRITE(*,*) "Before PrintGeom, L1059, Path.f90" |

1369 | 1 | pfleura2 | Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat) |

1370 | 1 | pfleura2 | STOP |

1371 | 1 | pfleura2 | END IF ! IF (Optgeom.GT.0) THEN |

1372 | 1 | pfleura2 | |

1373 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1374 | 1 | pfleura2 | ! End of GeomOpt |

1375 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1376 | 1 | pfleura2 | |

1377 | 1 | pfleura2 | IF (PathOnly) THEN |

1378 | 9 | pfleura2 | Call Header("PathOnly=.T. , Stop here") |

1379 | 1 | pfleura2 | Call Write_path(-1) |

1380 | 1 | pfleura2 | if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar) |

1381 | 1 | pfleura2 | STOP |

1382 | 1 | pfleura2 | END IF |

1383 | 1 | pfleura2 | |

1384 | 1 | pfleura2 | if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) |

1385 | 1 | pfleura2 | |

1386 | 1 | pfleura2 | DEALLOCATE(XyzGeomI,IntCoordI) |

1387 | 1 | pfleura2 | NGeomI=NGeomF |

1388 | 1 | pfleura2 | ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord)) |

1389 | 1 | pfleura2 | |

1390 | 1 | pfleura2 | IF (DEBUG) WRITE(*,*) ' Back from PathCreate, L965, Path.f90' |

1391 | 1 | pfleura2 | ! we print this path, which is the first one :-) called Iteration0 |

1392 | 1 | pfleura2 | ! we now have to calculate the gradient for each point (and the energy |

1393 | 1 | pfleura2 | ! if we work at 0K) |

1394 | 1 | pfleura2 | |

1395 | 1 | pfleura2 | if (DEBUG.AND.(COORD.EQ."BAKER")) THEN |

1396 | 1 | pfleura2 | DO IGeom=1,NGeomF |

1397 | 1 | pfleura2 | WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom |

1398 | 1 | pfleura2 | DO I=1,3*Nat-6 |

1399 | 1 | pfleura2 | WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) |

1400 | 1 | pfleura2 | END DO |

1401 | 1 | pfleura2 | END DO |

1402 | 1 | pfleura2 | END IF |

1403 | 1 | pfleura2 | |

1404 | 1 | pfleura2 | ! To calculate B^T^-1 for all extrapolated final geometries: |

1405 | 1 | pfleura2 | ! PFL 18 Jan 2008 -> |

1406 | 1 | pfleura2 | ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF |

1407 | 1 | pfleura2 | IF (COORD.EQ."BAKER") Call Calc_baker_allGeomF() |

1408 | 1 | pfleura2 | ! <-- PFL 18 Jan 2008 |

1409 | 1 | pfleura2 | |

1410 | 1 | pfleura2 | if (DEBUG.AND.(COORD.EQ."BAKER")) THEN |

1411 | 1 | pfleura2 | DO IGeom=1,NGeomF |

1412 | 1 | pfleura2 | WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom |

1413 | 1 | pfleura2 | DO I=1,3*Nat-6 |

1414 | 1 | pfleura2 | WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) |

1415 | 1 | pfleura2 | END DO |

1416 | 1 | pfleura2 | END DO |

1417 | 1 | pfleura2 | END IF |

1418 | 1 | pfleura2 | |

1419 | 1 | pfleura2 | IOpt=0 |

1420 | 2 | pfleura2 | Call EgradPath(IOpt) |

1421 | 1 | pfleura2 | |

1422 | 1 | pfleura2 | if (debug) WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt |

1423 | 1 | pfleura2 | Call Write_path(IOpt) |

1424 | 7 | pfleura2 | ! Shall we analyze the geometries ? |

1425 | 7 | pfleura2 | IF (AnaGeom) Call AnaPath(Iopt,IoDat) |

1426 | 7 | pfleura2 | |

1427 | 1 | pfleura2 | if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) |

1428 | 1 | pfleura2 | |

1429 | 1 | pfleura2 | ! We have the geometries, the first gradients... we will generate the first hessian matrices :-) |

1430 | 1 | pfleura2 | ! ... v3.71 Only if IniHUp=.TRUE. which is not the default. In this version, we take the hessian |

1431 | 1 | pfleura2 | ! of end points as Unity matrices, and we update them along the path. |

1432 | 1 | pfleura2 | |

1433 | 1 | pfleura2 | ! By default, Hess=Id |

1434 | 1 | pfleura2 | Hess=0. |

1435 | 1 | pfleura2 | IF(Coord .EQ. "ZMAT") Then |

1436 | 1 | pfleura2 | ! We use the fact that we know that approximate good hessian values |

1437 | 1 | pfleura2 | ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles |

1438 | 1 | pfleura2 | DO IGeom=1,NGeomF |

1439 | 1 | pfleura2 | Hess(IGeom,1,1)=0.5d0 |

1440 | 1 | pfleura2 | Hess(IGeom,2,2)=0.5d0 |

1441 | 1 | pfleura2 | Hess(IGeom,3,3)=0.1d0 |

1442 | 1 | pfleura2 | DO J=1,Nat-3 |

1443 | 1 | pfleura2 | Hess(IGeom,3*J+1,3*J+1)=0.5d0 |

1444 | 1 | pfleura2 | Hess(IGeom,3*J+2,3*J+2)=0.1d0 |

1445 | 1 | pfleura2 | Hess(IGeom,3*J+3,3*J+3)=0.02d0 |

1446 | 1 | pfleura2 | END DO |

1447 | 1 | pfleura2 | IF (HInv) THEN |

1448 | 1 | pfleura2 | DO I=1,NCoord |

1449 | 1 | pfleura2 | Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I) |

1450 | 1 | pfleura2 | END DO |

1451 | 1 | pfleura2 | END IF |

1452 | 1 | pfleura2 | END DO |

1453 | 1 | pfleura2 | ELSE |

1454 | 1 | pfleura2 | DO I=1,NCoord |

1455 | 1 | pfleura2 | DO J=1,NGeomF |

1456 | 1 | pfleura2 | Hess(J,I,I)=1.d0 |

1457 | 1 | pfleura2 | END DO |

1458 | 1 | pfleura2 | END DO |

1459 | 1 | pfleura2 | END IF |

1460 | 1 | pfleura2 | |

1461 | 1 | pfleura2 | IF (COORD .EQ. "BAKER") THEN |

1462 | 1 | pfleura2 | ! UMat(NPrim,NCoord) |

1463 | 1 | pfleura2 | ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6)) |

1464 | 1 | pfleura2 | Hprim=0.d0 |

1465 | 1 | pfleura2 | ScanCoord=>Coordinate |

1466 | 1 | pfleura2 | I=0 |

1467 | 1 | pfleura2 | DO WHILE (Associated(ScanCoord%next)) |

1468 | 1 | pfleura2 | I=I+1 |

1469 | 1 | pfleura2 | SELECT CASE (ScanCoord%Type) |

1470 | 1 | pfleura2 | CASE ('BOND') |

1471 | 1 | pfleura2 | !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF?? |

1472 | 1 | pfleura2 | Hprim(I,I) = 0.5d0 |

1473 | 1 | pfleura2 | !END DO |

1474 | 1 | pfleura2 | CASE ('ANGLE') |

1475 | 1 | pfleura2 | Hprim(I,I) = 0.2d0 |

1476 | 1 | pfleura2 | CASE ('DIHEDRAL') |

1477 | 1 | pfleura2 | Hprim(I,I) = 0.1d0 |

1478 | 1 | pfleura2 | END SELECT |

1479 | 1 | pfleura2 | ScanCoord => ScanCoord%next |

1480 | 1 | pfleura2 | END DO |

1481 | 1 | pfleura2 | |

1482 | 1 | pfleura2 | ! HprimU = H^prim U: |

1483 | 1 | pfleura2 | HprimU=0.d0 |

1484 | 1 | pfleura2 | DO I=1, 3*Nat-6 |

1485 | 1 | pfleura2 | DO J=1,NPrim |

1486 | 1 | pfleura2 | HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I) |

1487 | 1 | pfleura2 | END DO |

1488 | 1 | pfleura2 | END DO |

1489 | 1 | pfleura2 | |

1490 | 1 | pfleura2 | ! Hess = U^T HprimU = U^T H^prim U: |

1491 | 1 | pfleura2 | Hess=0.d0 |

1492 | 1 | pfleura2 | DO I=1, 3*Nat-6 |

1493 | 1 | pfleura2 | DO J=1,NPrim |

1494 | 1 | pfleura2 | ! UMat^T is needed below. |

1495 | 1 | pfleura2 | Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I) |

1496 | 1 | pfleura2 | END DO |

1497 | 1 | pfleura2 | END DO |

1498 | 1 | pfleura2 | DO K=2,NGeomF |

1499 | 1 | pfleura2 | Hess(K,:,:)=Hess(1,:,:) |

1500 | 1 | pfleura2 | END DO |

1501 | 1 | pfleura2 | DEALLOCATE(Hprim,HprimU) |

1502 | 1 | pfleura2 | end if ! ends IF COORD=BAKER |

1503 | 1 | pfleura2 | ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord)) |

1504 | 1 | pfleura2 | ALLOCATE(HessTmp(NCoord,NCoord)) |

1505 | 1 | pfleura2 | |

1506 | 1 | pfleura2 | IF (IniHUp) THEN |

1507 | 1 | pfleura2 | IF (FCalcHess) THEN |

1508 | 1 | pfleura2 | ! The numerical calculation of the Hessian has been switched on |

1509 | 1 | pfleura2 | DO IGeom=1,NGeomF |

1510 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1511 | 1 | pfleura2 | GeomTmp=IntCoordF(IGeom,:) |

1512 | 1 | pfleura2 | ELSE |

1513 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) |

1514 | 1 | pfleura2 | END IF |

1515 | 1 | pfleura2 | Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt) |

1516 | 1 | pfleura2 | END DO |

1517 | 1 | pfleura2 | ELSE |

1518 | 1 | pfleura2 | ! We generate a forward hessian and a backward one and we mix them. |

1519 | 1 | pfleura2 | ! First, the forward way: |

1520 | 1 | pfleura2 | DO IGeom=2,NGeomF-1 |

1521 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1522 | 1 | pfleura2 | GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:) |

1523 | 1 | pfleura2 | ELSE |

1524 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/)) |

1525 | 1 | pfleura2 | END IF |

1526 | 1 | pfleura2 | GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:) |

1527 | 1 | pfleura2 | HessTmp=Hess(IGeom-1,:,:) |

1528 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

1529 | 1 | pfleura2 | Hess(IGeom,:,:)=HessTmp |

1530 | 1 | pfleura2 | END DO |

1531 | 1 | pfleura2 | |

1532 | 1 | pfleura2 | ! Now, the backward way: |

1533 | 1 | pfleura2 | HessTmp=Hess(NGeomF,:,:) |

1534 | 1 | pfleura2 | DO IGeom=NGeomF-1,2,-1 |

1535 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1536 | 1 | pfleura2 | GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:) |

1537 | 1 | pfleura2 | ELSE |

1538 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) |

1539 | 1 | pfleura2 | END IF |

1540 | 1 | pfleura2 | GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:) |

1541 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

1542 | 1 | pfleura2 | alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.) |

1543 | 1 | pfleura2 | ca=cos(alpha) |

1544 | 1 | pfleura2 | sa=sin(alpha) |

1545 | 1 | pfleura2 | Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:) |

1546 | 1 | pfleura2 | END DO |

1547 | 1 | pfleura2 | END IF ! matches IF (FCalcHess) |

1548 | 1 | pfleura2 | END IF ! matches IF (IniHUp) THEN |

1549 | 1 | pfleura2 | |

1550 | 1 | pfleura2 | ! For debug purposes, we diagonalize the Hessian matrices |

1551 | 1 | pfleura2 | IF (debug) THEN |

1552 | 1 | pfleura2 | !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices" |

1553 | 1 | pfleura2 | DO I=1,NGeomF |

1554 | 1 | pfleura2 | WRITE(*,*) "DBG Main - Point ",I |

1555 | 1 | pfleura2 | Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord) |

1556 | 1 | pfleura2 | DO J=1,NCoord |

1557 | 1 | pfleura2 | WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord)) |

1558 | 1 | pfleura2 | END DO |

1559 | 1 | pfleura2 | END DO |

1560 | 1 | pfleura2 | END IF ! matches IF (debug) THEN |

1561 | 1 | pfleura2 | |

1562 | 1 | pfleura2 | ! we have the hessian matrices and the gradients, we can start the optimizations ! |

1563 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1564 | 1 | pfleura2 | ! |

1565 | 1 | pfleura2 | ! Main LOOP : Optimization of the Path |

1566 | 1 | pfleura2 | ! |

1567 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1568 | 1 | pfleura2 | if (debug) Print *, 'NGeomF=', NGeomF |

1569 | 1 | pfleura2 | allocation_flag = .TRUE. |

1570 | 1 | pfleura2 | |

1571 | 1 | pfleura2 | Fini=.FALSE. |

1572 | 1 | pfleura2 | |

1573 | 1 | pfleura2 | DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini)) |

1574 | 1 | pfleura2 | if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt |

1575 | 1 | pfleura2 | IOpt=IOpt+1 |

1576 | 1 | pfleura2 | |

1577 | 1 | pfleura2 | SELECT CASE (COORD) |

1578 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1579 | 1 | pfleura2 | GeomOld_all=IntCoordF |

1580 | 1 | pfleura2 | CASE ('CART','HYBRID') |

1581 | 1 | pfleura2 | GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/)) |

1582 | 1 | pfleura2 | CASE DEFAULT |

1583 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP' |

1584 | 1 | pfleura2 | STOP |

1585 | 1 | pfleura2 | END SELECT |

1586 | 1 | pfleura2 | |

1587 | 1 | pfleura2 | IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. & |

1588 | 1 | pfleura2 | valid("printtangent")) THEN |

1589 | 1 | pfleura2 | WRITE(*,*) "Visualization of tangents" |

1590 | 1 | pfleura2 | Idx=min(12,NCoord) |

1591 | 1 | pfleura2 | DO I=1,NGeomF |

1592 | 1 | pfleura2 | WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx) |

1593 | 1 | pfleura2 | WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx) |

1594 | 1 | pfleura2 | WRITE(*,*) |

1595 | 1 | pfleura2 | END DO |

1596 | 1 | pfleura2 | WRITE(*,*) "END of tangents" |

1597 | 1 | pfleura2 | END IF |

1598 | 1 | pfleura2 | |

1599 | 1 | pfleura2 | IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. & |

1600 | 1 | pfleura2 | valid("printgrad")) THEN |

1601 | 1 | pfleura2 | WRITE(*,*) "Visualization of gradients" |

1602 | 1 | pfleura2 | Idx=min(12,NCoord) |

1603 | 1 | pfleura2 | DO I=1,NGeomF |

1604 | 1 | pfleura2 | WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx) |

1605 | 1 | pfleura2 | WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx) |

1606 | 1 | pfleura2 | WRITe(*,*) |

1607 | 1 | pfleura2 | END DO |

1608 | 1 | pfleura2 | WRITE(*,*) "END of gradients" |

1609 | 1 | pfleura2 | END IF |

1610 | 1 | pfleura2 | |

1611 | 1 | pfleura2 | Fini=.TRUE. |

1612 | 1 | pfleura2 | IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE. |

1613 | 1 | pfleura2 | IGeom=1 |

1614 | 1 | pfleura2 | if (debug) WRITE(*,*) '**************************************' |

1615 | 1 | pfleura2 | if (debug) WRITE(*,*) 'Optimizing reactant' |

1616 | 1 | pfleura2 | if (debug) WRITE(*,*) '**************************************' |

1617 | 1 | pfleura2 | SELECT CASE (COORD) |

1618 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1619 | 1 | pfleura2 | SELECT CASE (Step_method) |

1620 | 1 | pfleura2 | CASE ('RFO') |

1621 | 1 | pfleura2 | !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT). |

1622 | 1 | pfleura2 | Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec) |

1623 | 1 | pfleura2 | Print *, GradTmp |

1624 | 1 | pfleura2 | CASE ('GDIIS') |

1625 | 1 | pfleura2 | ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).: |

1626 | 1 | pfleura2 | Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,& |

1627 | 1 | pfleura2 | Hess(1,:,:),NCoord,allocation_flag,ZeroVec) |

1628 | 1 | pfleura2 | Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)=' |

1629 | 1 | pfleura2 | Print *, GradTmp |

1630 | 1 | pfleura2 | CASE ('GEDIIS') |

1631 | 1 | pfleura2 | ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).: |

1632 | 1 | pfleura2 | ! Energies are updated in EgradPath.f90 |

1633 | 1 | pfleura2 | Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),& |

1634 | 1 | pfleura2 | NCoord,allocation_flag,ZeroVec) |

1635 | 1 | pfleura2 | !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag) |

1636 | 1 | pfleura2 | !GradTmp = GeomTmp - IntCoordF(1,:) |

1637 | 1 | pfleura2 | !Print *, 'Old Geometry:' |

1638 | 1 | pfleura2 | !Print *, IntCoordF(1,:) |

1639 | 1 | pfleura2 | Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)=' |

1640 | 1 | pfleura2 | Print *, GradTmp |

1641 | 1 | pfleura2 | !Print *, 'GeomTmp:' |

1642 | 1 | pfleura2 | !Print *, GeomTmp |

1643 | 1 | pfleura2 | GeomTmp(:) = GradTmp(:) + IntCoordF(1,:) |

1644 | 1 | pfleura2 | CASE DEFAULT |

1645 | 1 | pfleura2 | WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop" |

1646 | 1 | pfleura2 | STOP |

1647 | 1 | pfleura2 | END SELECT |

1648 | 1 | pfleura2 | CASE ('HYBRID','CART') |

1649 | 1 | pfleura2 | Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec) |

1650 | 1 | pfleura2 | CASE DEFAULT |

1651 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP' |

1652 | 1 | pfleura2 | STOP |

1653 | 1 | pfleura2 | END SELECT |

1654 | 1 | pfleura2 | |

1655 | 1 | pfleura2 | IF (debug) THEN |

1656 | 1 | pfleura2 | WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1' |

1657 | 1 | pfleura2 | Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat) |

1658 | 1 | pfleura2 | END IF |

1659 | 1 | pfleura2 | |

1660 | 1 | pfleura2 | ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax |

1661 | 1 | pfleura2 | NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) |

1662 | 1 | pfleura2 | FactStep=min(1.d0,MaxStep(Igeom)/NormStep) |

1663 | 1 | pfleura2 | Fini=(Fini.AND.(NormStep.LE.SThresh)) |

1664 | 1 | pfleura2 | OptReac=(NormStep.GT.SThresh) |

1665 | 1 | pfleura2 | IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac |

1666 | 1 | pfleura2 | NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/)))) |

1667 | 1 | pfleura2 | Fini=(Fini.AND.(NormGrad.LE.GThresh)) |

1668 | 1 | pfleura2 | OptReac=(OptReac.OR.(NormGrad.GT.GThresh)) |

1669 | 1 | pfleura2 | !Print *, 'Grad(1,:):' |

1670 | 1 | pfleura2 | !Print *, Grad(1,:) |

1671 | 1 | pfleura2 | IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac |

1672 | 1 | pfleura2 | |

1673 | 1 | pfleura2 | GradTmp=GradTmp*FactStep |

1674 | 1 | pfleura2 | |

1675 | 1 | pfleura2 | ! we take care that frozen atoms do not move. |

1676 | 1 | pfleura2 | IF (NFroz.NE.0) THEN |

1677 | 1 | pfleura2 | SELECT CASE (COORD) |

1678 | 1 | pfleura2 | CASE ('ZMAT','MIXED') |

1679 | 1 | pfleura2 | GradTmp(1:IntFroz)=0.d0 |

1680 | 1 | pfleura2 | CASE ('CART','HYBRID') |

1681 | 1 | pfleura2 | DO I=1,Nat |

1682 | 1 | pfleura2 | IF (FrozAtoms(I)) THEN |

1683 | 1 | pfleura2 | Iat=Order(I) |

1684 | 1 | pfleura2 | GradTmp(3*Iat-2:3*Iat)=0.d0 |

1685 | 1 | pfleura2 | END IF |

1686 | 1 | pfleura2 | END DO |

1687 | 1 | pfleura2 | CASE('BAKER') |

1688 | 1 | pfleura2 | GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..) |

1689 | 1 | pfleura2 | CASE DEFAULT |

1690 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP' |

1691 | 1 | pfleura2 | STOP |

1692 | 1 | pfleura2 | END SELECT |

1693 | 1 | pfleura2 | END IF ! matches IF (NFroz.NE.0) THEN |

1694 | 1 | pfleura2 | |

1695 | 1 | pfleura2 | IF (debug) THEN |

1696 | 1 | pfleura2 | !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1' |

1697 | 1 | pfleura2 | !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat) |

1698 | 1 | pfleura2 | END IF |

1699 | 1 | pfleura2 | |

1700 | 1 | pfleura2 | IF (debug) THEN |

1701 | 1 | pfleura2 | !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1' |

1702 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1703 | 1 | pfleura2 | WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree)) |

1704 | 1 | pfleura2 | ELSE |

1705 | 1 | pfleura2 | DO Iat=1,Nat |

1706 | 1 | pfleura2 | WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & |

1707 | 1 | pfleura2 | XyzGeomF(IGeom,1:3,Iat) |

1708 | 1 | pfleura2 | END DO |

1709 | 1 | pfleura2 | END IF |

1710 | 1 | pfleura2 | END IF ! matches IF (debug) THEN |

1711 | 1 | pfleura2 | |

1712 | 1 | pfleura2 | SELECT CASE (COORD) |

1713 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1714 | 1 | pfleura2 | IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1 |

1715 | 1 | pfleura2 | CASE ('HYBRID','CART') |

1716 | 1 | pfleura2 | XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/)) |

1717 | 1 | pfleura2 | CASE DEFAULT |

1718 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP' |

1719 | 1 | pfleura2 | STOP |

1720 | 1 | pfleura2 | END SELECT |

1721 | 1 | pfleura2 | |

1722 | 1 | pfleura2 | IF (debug) THEN |

1723 | 1 | pfleura2 | WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1' |

1724 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1725 | 1 | pfleura2 | WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree)) |

1726 | 1 | pfleura2 | ELSE |

1727 | 1 | pfleura2 | DO Iat=1,Nat |

1728 | 1 | pfleura2 | |

1729 | 1 | pfleura2 | !!!!!!!!!! There was an OrderInv here... should I put it back ? |

1730 | 1 | pfleura2 | ! It was with indzmat. IndZmat cannot appear here... |

1731 | 1 | pfleura2 | WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & |

1732 | 1 | pfleura2 | XyzGeomF(IGeom,1:3,Iat) |

1733 | 1 | pfleura2 | END DO |

1734 | 1 | pfleura2 | END IF |

1735 | 1 | pfleura2 | END IF ! matches IF (debug) THEN |

1736 | 1 | pfleura2 | |

1737 | 1 | pfleura2 | IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized" |

1738 | 1 | pfleura2 | ELSE ! Optreac |

1739 | 1 | pfleura2 | IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini |

1740 | 1 | pfleura2 | END IF ! matches IF (OptReac) THEN |

1741 | 1 | pfleura2 | |

1742 | 1 | pfleura2 | IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE. |

1743 | 1 | pfleura2 | IGeom=NGeomF |

1744 | 1 | pfleura2 | if (debug) WRITE(*,*) '******************' |

1745 | 1 | pfleura2 | if (debug) WRITE(*,*) 'Optimizing product' |

1746 | 1 | pfleura2 | if (debug) WRITE(*,*) '******************' |

1747 | 1 | pfleura2 | SELECT CASE (COORD) |

1748 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1749 | 1 | pfleura2 | Print *, 'Optimizing product' |

1750 | 1 | pfleura2 | SELECT CASE (Step_method) |

1751 | 1 | pfleura2 | CASE ('RFO') |

1752 | 1 | pfleura2 | !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT). |

1753 | 1 | pfleura2 | Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec) |

1754 | 1 | pfleura2 | Print *, GradTmp |

1755 | 1 | pfleura2 | CASE ('GDIIS') |

1756 | 1 | pfleura2 | ! GradTmp becomes "step" below and has INTENT(OUT): |

1757 | 1 | pfleura2 | Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),& |

1758 | 1 | pfleura2 | HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec) |

1759 | 1 | pfleura2 | Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)=' |

1760 | 1 | pfleura2 | Print *, GradTmp |

1761 | 1 | pfleura2 | CASE ('GEDIIS') |

1762 | 1 | pfleura2 | ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT): |

1763 | 1 | pfleura2 | Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),& |

1764 | 1 | pfleura2 | Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec) |

1765 | 1 | pfleura2 | Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)=' |

1766 | 1 | pfleura2 | Print *, GradTmp |

1767 | 1 | pfleura2 | CASE DEFAULT |

1768 | 1 | pfleura2 | WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop" |

1769 | 1 | pfleura2 | STOP |

1770 | 1 | pfleura2 | END SELECT |

1771 | 1 | pfleura2 | CASE ('HYBRID','CART') |

1772 | 1 | pfleura2 | Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec) |

1773 | 1 | pfleura2 | CASE DEFAULT |

1774 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP' |

1775 | 1 | pfleura2 | STOP |

1776 | 1 | pfleura2 | END SELECT |

1777 | 1 | pfleura2 | |

1778 | 1 | pfleura2 | ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax |

1779 | 1 | pfleura2 | NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) |

1780 | 1 | pfleura2 | FactStep=min(1.d0,MaxStep(IGeom)/NormStep) |

1781 | 1 | pfleura2 | Fini=(Fini.AND.(NormStep.LE.SThresh)) |

1782 | 1 | pfleura2 | OptProd=.NOT.(NormStep.LE.SThresh) |

1783 | 1 | pfleura2 | NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/)))) |

1784 | 1 | pfleura2 | Fini=(Fini.AND.(NormGrad.LE.GThresh)) |

1785 | 1 | pfleura2 | OptProd=(OptProd.OR.(NormGrad.GT.GThresh)) |

1786 | 1 | pfleura2 | IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd |

1787 | 1 | pfleura2 | |

1788 | 1 | pfleura2 | GradTmp=GradTmp*FactStep |

1789 | 1 | pfleura2 | |

1790 | 1 | pfleura2 | ! we take care that frozen atoms do not move |

1791 | 1 | pfleura2 | IF (NFroz.NE.0) THEN |

1792 | 1 | pfleura2 | SELECT CASE (COORD) |

1793 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1794 | 1 | pfleura2 | GradTmp(1:IntFroz)=0.d0 |

1795 | 1 | pfleura2 | CASE ('CART','HYBRID') |

1796 | 1 | pfleura2 | DO I=1,Nat |

1797 | 1 | pfleura2 | IF (FrozAtoms(I)) THEN |

1798 | 1 | pfleura2 | Iat=Order(I) |

1799 | 1 | pfleura2 | GradTmp(3*Iat-2:3*Iat)=0.d0 |

1800 | 1 | pfleura2 | END IF |

1801 | 1 | pfleura2 | END DO |

1802 | 1 | pfleura2 | CASE DEFAULT |

1803 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP' |

1804 | 1 | pfleura2 | STOP |

1805 | 1 | pfleura2 | END SELECT |

1806 | 1 | pfleura2 | END IF ! matches IF (NFroz.NE.0) THEN |

1807 | 1 | pfleura2 | |

1808 | 1 | pfleura2 | IF (debug) THEN |

1809 | 1 | pfleura2 | WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt |

1810 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1811 | 1 | pfleura2 | WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) |

1812 | 1 | pfleura2 | ELSE |

1813 | 1 | pfleura2 | DO Iat=1,Nat |

1814 | 1 | pfleura2 | WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & |

1815 | 1 | pfleura2 | XyzGeomF(IGeom,1:3,Iat) |

1816 | 1 | pfleura2 | END DO |

1817 | 1 | pfleura2 | END IF |

1818 | 1 | pfleura2 | END IF |

1819 | 1 | pfleura2 | |

1820 | 1 | pfleura2 | SELECT CASE (COORD) |

1821 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1822 | 1 | pfleura2 | IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF. |

1823 | 1 | pfleura2 | CASE ('HYBRID','CART') |

1824 | 1 | pfleura2 | XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) |

1825 | 1 | pfleura2 | CASE DEFAULT |

1826 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP' |

1827 | 1 | pfleura2 | STOP |

1828 | 1 | pfleura2 | END SELECT |

1829 | 1 | pfleura2 | |

1830 | 1 | pfleura2 | IF (debug) THEN |

1831 | 1 | pfleura2 | WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt |

1832 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1833 | 1 | pfleura2 | WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) |

1834 | 1 | pfleura2 | ELSE |

1835 | 1 | pfleura2 | DO Iat=1,Nat |

1836 | 1 | pfleura2 | WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & |

1837 | 1 | pfleura2 | XyzGeomF(IGeom,1:3,Iat) |

1838 | 1 | pfleura2 | END DO |

1839 | 1 | pfleura2 | END IF |

1840 | 1 | pfleura2 | END IF |

1841 | 1 | pfleura2 | ELSE ! Optprod |

1842 | 1 | pfleura2 | IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini |

1843 | 1 | pfleura2 | END IF !matches IF (OptProd) THEN |

1844 | 1 | pfleura2 | |

1845 | 1 | pfleura2 | IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini |

1846 | 1 | pfleura2 | |

1847 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1848 | 1 | pfleura2 | ! |

1849 | 1 | pfleura2 | ! Optimizing other geometries |

1850 | 1 | pfleura2 | ! |

1851 | 1 | pfleura2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

1852 | 1 | pfleura2 | |

1853 | 1 | pfleura2 | |

1854 | 1 | pfleura2 | |

1855 | 1 | pfleura2 | DO IGeom=2,NGeomF-1 ! matches at L1556 |

1856 | 1 | pfleura2 | if (debug) WRITE(*,*) '****************************' |

1857 | 1 | pfleura2 | if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom |

1858 | 1 | pfleura2 | if (debug) WRITE(*,*) '****************************' |

1859 | 1 | pfleura2 | |

1860 | 1 | pfleura2 | SELECT CASE (COORD) |

1861 | 1 | pfleura2 | CASE ('ZMAT','BAKER','MIXED') |

1862 | 1 | pfleura2 | GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/)) |

1863 | 1 | pfleura2 | ! PFL 6 Apr 2008 -> |

1864 | 1 | pfleura2 | ! Special case, if FTan=0. we do not consider Tangent at all. |

1865 | 1 | pfleura2 | ! To DO: add the full treatment of FTan |

1866 | 1 | pfleura2 | if (FTan==0.) GradTmp2=ZeroVec |

1867 | 1 | pfleura2 | ! <- PFL 6 Apr 2008 |

1868 | 1 | pfleura2 | if (debug) THEN |

1869 | 1 | pfleura2 | Print *, 'L1500, IntTangent(',IGeom,',:)=' |

1870 | 1 | pfleura2 | Print *, IntTangent(IGeom,:) |

1871 | 1 | pfleura2 | END IF |

1872 | 1 | pfleura2 | !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT). |

1873 | 1 | pfleura2 | SELECT CASE (Step_method) |

1874 | 1 | pfleura2 | CASE ('RFO') |

1875 | 1 | pfleura2 | Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2) |

1876 | 1 | pfleura2 | CASE ('GDIIS') |

1877 | 1 | pfleura2 | !The following has no effect at IOpt==1 |

1878 | 1 | pfleura2 | !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom |

1879 | 1 | pfleura2 | Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,& |

1880 | 1 | pfleura2 | Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2) |

1881 | 1 | pfleura2 | Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)=' |

1882 | 1 | pfleura2 | Print *, GradTmp |

1883 | 1 | pfleura2 | CASE ('GEDIIS') |

1884 | 1 | pfleura2 | ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT): |

1885 | 1 | pfleura2 | Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),& |

1886 | 1 | pfleura2 | Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2) |

1887 | 1 | pfleura2 | Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)=' |

1888 | 1 | pfleura2 | Print *, GradTmp |

1889 | 1 | pfleura2 | CASE DEFAULT |

1890 | 1 | pfleura2 | WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop" |

1891 | 1 | pfleura2 | STOP |

1892 | 1 | pfleura2 | END SELECT |

1893 | 1 | pfleura2 | CASE ('HYBRID','CART') |

1894 | 1 | pfleura2 | ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat |

1895 | 1 | pfleura2 | ! so we change it |

1896 | 1 | pfleura2 | DO I=1,Nat |

1897 | 1 | pfleura2 | DO J=1,3 |

1898 | 1 | pfleura2 | GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I) |

1899 | 1 | pfleura2 | END DO |

1900 | 1 | pfleura2 | END DO |

1901 | 1 | pfleura2 | ! GradTmp2=XyzTangent(IGeom,:) |

1902 | 1 | pfleura2 | ! PFL 6 Apr 2008 -> |

1903 | 1 | pfleura2 | ! Special case, if FTan=0. we do not consider Tangent at all. |

1904 | 1 | pfleura2 | ! To DO: add the full treatment of FTan |

1905 | 1 | pfleura2 | if (FTan==0.) GradTmp2=ZeroVec |

1906 | 1 | pfleura2 | ! <- PFL 6 Apr 2008 |

1907 | 1 | pfleura2 | |

1908 | 1 | pfleura2 | Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2) |

1909 | 1 | pfleura2 | CASE DEFAULT |

1910 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP' |

1911 | 1 | pfleura2 | STOP |

1912 | 1 | pfleura2 | END SELECT |

1913 | 1 | pfleura2 | |

1914 | 1 | pfleura2 | ! we take care that frozen atoms do not move |

1915 | 1 | pfleura2 | IF (NFroz.NE.0) THEN |

1916 | 1 | pfleura2 | SELECT CASE (COORD) |

1917 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1918 | 1 | pfleura2 | IF (debug) THEN |

1919 | 1 | pfleura2 | WRITE(*,*) 'Step computed. Coord=',Coord |

1920 | 1 | pfleura2 | WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz |

1921 | 1 | pfleura2 | END IF |

1922 | 1 | pfleura2 | GradTmp(1:IntFroz)=0.d0 |

1923 | 1 | pfleura2 | CASE ('CART','HYBRID') |

1924 | 1 | pfleura2 | if (debug) THEN |

1925 | 1 | pfleura2 | WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing" |

1926 | 1 | pfleura2 | DO I=1,Nat |

1927 | 1 | pfleura2 | WRITe(*,*) I,GradTmp(3*I-2:3*I) |

1928 | 1 | pfleura2 | END DO |

1929 | 1 | pfleura2 | END IF |

1930 | 1 | pfleura2 | DO I=1,Nat |

1931 | 1 | pfleura2 | IF (FrozAtoms(I)) THEN |

1932 | 1 | pfleura2 | if (debug) THEN |

1933 | 1 | pfleura2 | write(*,*) 'Step Computed. Coord=',Coord |

1934 | 1 | pfleura2 | WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I) |

1935 | 1 | pfleura2 | END IF |

1936 | 1 | pfleura2 | Iat=Order(I) |

1937 | 1 | pfleura2 | GradTmp(3*Iat-2:3*Iat)=0.d0 |

1938 | 1 | pfleura2 | END IF |

1939 | 1 | pfleura2 | END DO |

1940 | 1 | pfleura2 | |

1941 | 1 | pfleura2 | if (debug) THEN |

1942 | 1 | pfleura2 | WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing" |

1943 | 1 | pfleura2 | DO I=1,Nat |

1944 | 1 | pfleura2 | WRITe(*,*) I,GradTmp(3*I-2:3*I) |

1945 | 1 | pfleura2 | END DO |

1946 | 1 | pfleura2 | END IF |

1947 | 1 | pfleura2 | |

1948 | 1 | pfleura2 | CASE DEFAULT |

1949 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP' |

1950 | 1 | pfleura2 | STOP |

1951 | 1 | pfleura2 | END SELECT |

1952 | 1 | pfleura2 | END IF !matches IF (NFroz.NE.0) THEN |

1953 | 1 | pfleura2 | |

1954 | 1 | pfleura2 | IF (debug) THEN |

1955 | 1 | pfleura2 | SELECT CASE (COORD) |

1956 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

1957 | 1 | pfleura2 | WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree)) |

1958 | 1 | pfleura2 | WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree)) |

1959 | 1 | pfleura2 | CASE ('CART','HYBRID') |

1960 | 1 | pfleura2 | WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree)) |

1961 | 1 | pfleura2 | WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree)) |

1962 | 1 | pfleura2 | CASE DEFAULT |

1963 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP' |

1964 | 1 | pfleura2 | STOP |

1965 | 1 | pfleura2 | END SELECT |

1966 | 1 | pfleura2 | END IF ! matches if (debug) THEN |

1967 | 1 | pfleura2 | |

1968 | 1 | pfleura2 | ! we check that the maximum displacement is not greater than Smax |

1969 | 1 | pfleura2 | If (debug) WRITE(*,*) "Fini before test:",Fini |

1970 | 1 | pfleura2 | NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) |

1971 | 1 | pfleura2 | FactStep=min(1.d0,maxStep(Igeom)/NormStep) |

1972 | 1 | pfleura2 | Fini=(Fini.AND.(NormStep.LE.SThresh)) |

1973 | 1 | pfleura2 | IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini |

1974 | 1 | pfleura2 | |

1975 | 1 | pfleura2 | GradTmp=GraDTmp*FactStep |

1976 | 1 | pfleura2 | |

1977 | 1 | pfleura2 | if (debug) WRITE(*,*) "DBG Main, L1475, FactStep=",FactStep |

1978 | 1 | pfleura2 | if (debug.AND.(COORD.EQ.'ZMAT')) WRITE(*,'(A,12(1X,F10.6))') 'Scaled Step:',GradTmp(1:min(12,NFree)) |

1979 | 1 | pfleura2 | |

1980 | 1 | pfleura2 | IF (debug) THEN |

1981 | 1 | pfleura2 | WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt |

1982 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

1983 | 1 | pfleura2 | WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) |

1984 | 1 | pfleura2 | ELSE |

1985 | 1 | pfleura2 | DO Iat=1,Nat |

1986 | 1 | pfleura2 | WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & |

1987 | 1 | pfleura2 | XyzGeomF(IGeom,1:3,Iat) |

1988 | 1 | pfleura2 | END DO |

1989 | 1 | pfleura2 | END IF |

1990 | 1 | pfleura2 | END IF ! matches if (debug) THEN |

1991 | 1 | pfleura2 | |

1992 | 1 | pfleura2 | SELECT CASE (COORD) |

1993 | 1 | pfleura2 | CASE ('ZMAT') |

1994 | 1 | pfleura2 | IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) |

1995 | 1 | pfleura2 | if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) |

1996 | 1 | pfleura2 | WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt |

1997 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) |

1998 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & |

1999 | 1 | pfleura2 | OrderInv(indzmat(2,2)),IntcoordF(IGeom,1) |

2000 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |

2001 | 1 | pfleura2 | OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi |

2002 | 1 | pfleura2 | Idx=4 |

2003 | 1 | pfleura2 | DO Iat=4,Nat |

2004 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |

2005 | 1 | pfleura2 | OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & |

2006 | 1 | pfleura2 | OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & |

2007 | 1 | pfleura2 | OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi |

2008 | 1 | pfleura2 | Idx=Idx+3 |

2009 | 1 | pfleura2 | END DO |

2010 | 1 | pfleura2 | |

2011 | 1 | pfleura2 | CASE ('BAKER') |

2012 | 1 | pfleura2 | IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) |

2013 | 1 | pfleura2 | if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) |

2014 | 1 | pfleura2 | WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt |

2015 | 1 | pfleura2 | WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull." |

2016 | 1 | pfleura2 | |

2017 | 1 | pfleura2 | CASE ('MIXED') |

2018 | 1 | pfleura2 | IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:) |

2019 | 1 | pfleura2 | if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) |

2020 | 1 | pfleura2 | WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt |

2021 | 1 | pfleura2 | DO Iat=1,NCart |

2022 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat) |

2023 | 1 | pfleura2 | END DO |

2024 | 1 | pfleura2 | Idx=3*NCart+1 |

2025 | 1 | pfleura2 | IF (NCart.EQ.1) THEN |

2026 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & |

2027 | 1 | pfleura2 | OrderInv(indzmat(2,2)),IntcoordF(IGeom,4) |

2028 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |

2029 | 1 | pfleura2 | OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi |

2030 | 1 | pfleura2 | |

2031 | 1 | pfleura2 | Idx=7 |

2032 | 1 | pfleura2 | END IF |

2033 | 1 | pfleura2 | IF (NCart.EQ.2) THEN |

2034 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |

2035 | 1 | pfleura2 | OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi |

2036 | 1 | pfleura2 | Idx=9 |

2037 | 1 | pfleura2 | END IF |

2038 | 1 | pfleura2 | |

2039 | 1 | pfleura2 | |

2040 | 1 | pfleura2 | DO Iat=max(NCart+1,4),Nat |

2041 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |

2042 | 1 | pfleura2 | OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), & |

2043 | 1 | pfleura2 | OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, & |

2044 | 1 | pfleura2 | OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi |

2045 | 1 | pfleura2 | Idx=Idx+3 |

2046 | 1 | pfleura2 | END DO |

2047 | 1 | pfleura2 | if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:) |

2048 | 1 | pfleura2 | |

2049 | 1 | pfleura2 | CASE ('HYBRID','CART') |

2050 | 1 | pfleura2 | XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/)) |

2051 | 1 | pfleura2 | CASE DEFAULT |

2052 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP' |

2053 | 1 | pfleura2 | STOP |

2054 | 1 | pfleura2 | END SELECT |

2055 | 1 | pfleura2 | |

2056 | 1 | pfleura2 | IF (debug) THEN |

2057 | 1 | pfleura2 | WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt |

2058 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

2059 | 1 | pfleura2 | WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree)) |

2060 | 1 | pfleura2 | ELSE |

2061 | 1 | pfleura2 | DO Iat=1,Nat |

2062 | 1 | pfleura2 | WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), & |

2063 | 1 | pfleura2 | XyzGeomF(IGeom,1:3,Iat) |

2064 | 1 | pfleura2 | END DO |

2065 | 1 | pfleura2 | END IF |

2066 | 1 | pfleura2 | END IF ! matches if (debug) THEN |

2067 | 1 | pfleura2 | |

2068 | 1 | pfleura2 | ! We project out the tangential part of the gradient to know if we are done |

2069 | 1 | pfleura2 | GradTmp=Grad(IGeom,:) |

2070 | 1 | pfleura2 | NormGrad=sqrt(dot_product(GradTmp,GradTmp)) |

2071 | 1 | pfleura2 | if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad |

2072 | 1 | pfleura2 | SELECT CASE (COORD) |

2073 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

2074 | 1 | pfleura2 | S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) |

2075 | 1 | pfleura2 | Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:)) |

2076 | 1 | pfleura2 | GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:) |

2077 | 1 | pfleura2 | Ca=S/(sqrt(Norm)*NormGrad) |

2078 | 1 | pfleura2 | S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:)) |

2079 | 1 | pfleura2 | GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:) |

2080 | 1 | pfleura2 | CASE ('CART','HYBRID') |

2081 | 1 | pfleura2 | S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) |

2082 | 1 | pfleura2 | Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:)) |

2083 | 1 | pfleura2 | Ca=S/(sqrt(Norm)*NormGrad) |

2084 | 1 | pfleura2 | GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:) |

2085 | 1 | pfleura2 | S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:)) |

2086 | 1 | pfleura2 | GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:) |

2087 | 1 | pfleura2 | CASE DEFAULT |

2088 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP' |

2089 | 1 | pfleura2 | STOP |

2090 | 1 | pfleura2 | END SELECT |

2091 | 1 | pfleura2 | IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi |

2092 | 1 | pfleura2 | NormGrad=sqrt(DOT_Product(GradTmp,GradTmp)) |

2093 | 1 | pfleura2 | Fini=(Fini.AND.(NormGrad.LE.GThresh)) |

2094 | 1 | pfleura2 | IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini |

2095 | 1 | pfleura2 | |

2096 | 1 | pfleura2 | END DO ! matches DO IGeom=2,NGeomF-1 |

2097 | 1 | pfleura2 | ! we save the old gradients |

2098 | 1 | pfleura2 | GradOld=Grad |

2099 | 1 | pfleura2 | EPathOld=Energies |

2100 | 1 | pfleura2 | |

2101 | 7 | pfleura2 | |

2102 | 7 | pfleura2 | |

2103 | 1 | pfleura2 | ! Shall we re-parameterize the path ? |

2104 | 1 | pfleura2 | ! PFL 2007/Apr/10 -> |

2105 | 1 | pfleura2 | ! We call PathCreate in all cases, it will take care of the |

2106 | 1 | pfleura2 | ! reparameterization, as well as calculating the tangents. |

2107 | 1 | pfleura2 | ! IF (MOD(IOpt,IReparam).EQ.0) THEN |

2108 | 1 | pfleura2 | XyzGeomI=XyzGeomF |

2109 | 1 | pfleura2 | IntCoordI=IntCoordF |

2110 | 1 | pfleura2 | |

2111 | 1 | pfleura2 | Call PathCreate(IOpt) |

2112 | 1 | pfleura2 | ! END IF |

2113 | 1 | pfleura2 | ! <- PFL 2007/Apr/10 |

2114 | 1 | pfleura2 | |

2115 | 1 | pfleura2 | if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)') 'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini |

2116 | 1 | pfleura2 | IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN |

2117 | 1 | pfleura2 | |

2118 | 1 | pfleura2 | ! We have the new path, we calculate its energy and gradient |

2119 | 1 | pfleura2 | |

2120 | 2 | pfleura2 | Call EgradPath(IOpt) |

2121 | 1 | pfleura2 | !IF(IOPT .EQ. 10) Then |

2122 | 1 | pfleura2 | ! Print *, 'Stopping at my checkpoint.' |

2123 | 1 | pfleura2 | ! STOP !This is my temporary checkpoint. |

2124 | 1 | pfleura2 | !ENDIF |

2125 | 1 | pfleura2 | |

2126 | 1 | pfleura2 | ! PFL 31 Mar 2008 -> |

2127 | 1 | pfleura2 | ! Taken from v3.99 but we added a flag... |

2128 | 1 | pfleura2 | ! Added in v3.99 : MaxStep is modified according to the evolution of energies |

2129 | 1 | pfleura2 | ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 |

2130 | 1 | pfleura2 | ! else it is decreased by 0.8 |

2131 | 1 | pfleura2 | |

2132 | 1 | pfleura2 | if (DynMaxStep) THEN |

2133 | 1 | pfleura2 | if (debug) WRITE(*,*) "Dynamically updating the Maximum step" |

2134 | 1 | pfleura2 | if (OptReac) THEN |

2135 | 1 | pfleura2 | If (Energies(1)<EPathOld(1)) Then |

2136 | 1 | pfleura2 | MaxStep(1)=min(1.1*MaxStep(1),2.*SMax) |

2137 | 1 | pfleura2 | ELSE |

2138 | 1 | pfleura2 | MaxStep(1)=max(0.8*MaxStep(1),SMax/2.) |

2139 | 1 | pfleura2 | END IF |

2140 | 1 | pfleura2 | END IF |

2141 | 1 | pfleura2 | |

2142 | 1 | pfleura2 | if (OptProd) THEN |

2143 | 1 | pfleura2 | If (Energies(NGeomF)<EPathOld(NGeomF)) Then |

2144 | 1 | pfleura2 | MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax) |

2145 | 1 | pfleura2 | ELSE |

2146 | 1 | pfleura2 | MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.) |

2147 | 1 | pfleura2 | END IF |

2148 | 1 | pfleura2 | END IF |

2149 | 1 | pfleura2 | |

2150 | 1 | pfleura2 | DO IGeom=2,NGeomF-1 |

2151 | 1 | pfleura2 | If (Energies(IGeom)<EPathOld(IGeom)) Then |

2152 | 1 | pfleura2 | MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax) |

2153 | 1 | pfleura2 | ELSE |

2154 | 1 | pfleura2 | MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.) |

2155 | 1 | pfleura2 | END IF |

2156 | 1 | pfleura2 | END DO |

2157 | 1 | pfleura2 | if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) |

2158 | 1 | pfleura2 | END IF ! test on DynMaxStep |

2159 | 1 | pfleura2 | if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF) |

2160 | 1 | pfleura2 | ! <- PFL 31 MAr 2008 |

2161 | 1 | pfleura2 | ! Also XyzGeomF is updated in EgradPath for Baker Case. |

2162 | 1 | pfleura2 | ! It should get updated for other cases too. |

2163 | 1 | pfleura2 | |

2164 | 1 | pfleura2 | DO IGeom=1,NGeomF |

2165 | 1 | pfleura2 | SELECT CASE (COORD) |

2166 | 1 | pfleura2 | CASE('ZMAT') |

2167 | 1 | pfleura2 | WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt |

2168 | 1 | pfleura2 | GeomTmp=IntCoordF(IGeom,:) |

2169 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1)))) |

2170 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & |

2171 | 1 | pfleura2 | OrderInv(indzmat(2,2)),Geomtmp(1) |

2172 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |

2173 | 1 | pfleura2 | OrderInv(indzmat(3,2)),Geomtmp(2), & |

2174 | 1 | pfleura2 | OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi |

2175 | 1 | pfleura2 | Idx=4 |

2176 | 1 | pfleura2 | DO Iat=4,Nat |

2177 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |

2178 | 1 | pfleura2 | OrderInv(indzmat(iat,2)),Geomtmp(Idx), & |

2179 | 1 | pfleura2 | OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, & |

2180 | 1 | pfleura2 | OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi |

2181 | 1 | pfleura2 | Idx=Idx+3 |

2182 | 1 | pfleura2 | END DO |

2183 | 1 | pfleura2 | CASE('BAKER') |

2184 | 1 | pfleura2 | !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt |

2185 | 1 | pfleura2 | !GeomTmp=IntCoordF(IGeom,:) |

2186 | 1 | pfleura2 | CASE('CART','HYBRID','MIXED') |

2187 | 1 | pfleura2 | WRITE(IOOUT,'(1X,I5)') Nat |

2188 | 1 | pfleura2 | WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt |

2189 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) |

2190 | 1 | pfleura2 | DO I=1,Nat |

2191 | 1 | pfleura2 | Iat=I |

2192 | 1 | pfleura2 | If (renum) Iat=OrderInv(I) |

2193 | 1 | pfleura2 | WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat) |

2194 | 1 | pfleura2 | END DO |

2195 | 1 | pfleura2 | CASE DEFAULT |

2196 | 1 | pfleura2 | WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP' |

2197 | 1 | pfleura2 | STOP |

2198 | 1 | pfleura2 | END SELECT |

2199 | 1 | pfleura2 | END DO ! matches DO IGeom=1,NGeomF |

2200 | 1 | pfleura2 | |

2201 | 1 | pfleura2 | Call Write_path(IOpt) |

2202 | 7 | pfleura2 | ! Shall we analyze the geometries ? |

2203 | 7 | pfleura2 | IF (AnaGeom) Call AnaPath(Iopt,IoDat) |

2204 | 1 | pfleura2 | |

2205 | 1 | pfleura2 | ! We update the Hessian matrices |

2206 | 1 | pfleura2 | IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord |

2207 | 1 | pfleura2 | ! First using the displacement from the old path |

2208 | 1 | pfleura2 | IGeom0=2 |

2209 | 1 | pfleura2 | IGeomF=NGeomF-1 |

2210 | 1 | pfleura2 | IF (OptReac) IGeom0=1 |

2211 | 1 | pfleura2 | IF (OptProd) IGeomF=NGeomF |

2212 | 1 | pfleura2 | |

2213 | 1 | pfleura2 | IF (FCalcHess) THEN |

2214 | 1 | pfleura2 | DO IGeom=IGeom0,IGeomF |

2215 | 1 | pfleura2 | SELECT CASE (COORD) |

2216 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

2217 | 1 | pfleura2 | GeomTmp2=IntCoordF(IGeom,:) |

2218 | 1 | pfleura2 | CASE ('CART','HYBRID') |

2219 | 1 | pfleura2 | GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) |

2220 | 1 | pfleura2 | CASE DEFAULT |

2221 | 1 | pfleura2 | WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" |

2222 | 1 | pfleura2 | STOP |

2223 | 1 | pfleura2 | END SELECT |

2224 | 1 | pfleura2 | Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt) |

2225 | 1 | pfleura2 | END DO |

2226 | 1 | pfleura2 | ELSE |

2227 | 1 | pfleura2 | DO IGeom=IGeom0,IGeomF |

2228 | 1 | pfleura2 | GeomTmp=GeomOld_all(IGeom,:) |

2229 | 1 | pfleura2 | SELECT CASE (COORD) |

2230 | 1 | pfleura2 | CASE ('ZMAT','MIXED','BAKER') |

2231 | 1 | pfleura2 | GeomTmp2=IntCoordF(IGeom,:) |

2232 | 1 | pfleura2 | CASE ('CART','HYBRID') |

2233 | 1 | pfleura2 | GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) |

2234 | 1 | pfleura2 | CASE DEFAULT |

2235 | 1 | pfleura2 | WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP" |

2236 | 1 | pfleura2 | STOP |

2237 | 1 | pfleura2 | END SELECT |

2238 | 1 | pfleura2 | GeomTmp=GeomTmp2-GeomTmp |

2239 | 1 | pfleura2 | GradTmp=Grad(IGeom,:)-GradOld(IGeom,:) |

2240 | 1 | pfleura2 | HessTmp=Hess(IGeom,:,:) |

2241 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

2242 | 1 | pfleura2 | Hess(IGeom,:,:)=HessTmp |

2243 | 1 | pfleura2 | END DO ! matches DO IGeom=IGeom0,IGeomF |

2244 | 1 | pfleura2 | |

2245 | 1 | pfleura2 | ! Second using the neighbour points |

2246 | 1 | pfleura2 | IF (HupNeighbour) THEN |

2247 | 1 | pfleura2 | ! Only one point for end points. |

2248 | 1 | pfleura2 | IF (OptReac) THEN |

2249 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

2250 | 1 | pfleura2 | GeomTmp=IntCoordF(1,:)-IntCoordF(2,:) |

2251 | 1 | pfleura2 | ELSE |

2252 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/)) |

2253 | 1 | pfleura2 | END IF |

2254 | 1 | pfleura2 | GradTmp=Grad(1,:)-Grad(2,:) |

2255 | 1 | pfleura2 | HessTmp=Hess(1,:,:) |

2256 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

2257 | 1 | pfleura2 | Hess(1,:,:)=HessTmp |

2258 | 1 | pfleura2 | END IF |

2259 | 1 | pfleura2 | |

2260 | 1 | pfleura2 | IF (OptProd) THEN |

2261 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

2262 | 1 | pfleura2 | GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:) |

2263 | 1 | pfleura2 | ELSE |

2264 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/)) |

2265 | 1 | pfleura2 | END IF |

2266 | 1 | pfleura2 | GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:) |

2267 | 1 | pfleura2 | HessTmp=Hess(NGeomF,:,:) |

2268 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

2269 | 1 | pfleura2 | Hess(NGeomF,:,:)=HessTmp |

2270 | 1 | pfleura2 | END IF |

2271 | 1 | pfleura2 | |

2272 | 1 | pfleura2 | ! Two points for the rest of the path |

2273 | 1 | pfleura2 | DO IGeom=2,NGeomF-1 |

2274 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN |

2275 | 1 | pfleura2 | GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:) |

2276 | 1 | pfleura2 | ELSE |

2277 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/)) |

2278 | 1 | pfleura2 | END IF |

2279 | 1 | pfleura2 | GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:) |

2280 | 1 | pfleura2 | HessTmp=Hess(IGeom,:,:) |

2281 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

2282 | 1 | pfleura2 | |

2283 | 1 | pfleura2 | IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN |

2284 | 1 | pfleura2 | GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:) |

2285 | 1 | pfleura2 | ELSE |

2286 | 1 | pfleura2 | GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/)) |

2287 | 1 | pfleura2 | END IF |

2288 | 1 | pfleura2 | GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:) |

2289 | 1 | pfleura2 | |

2290 | 1 | pfleura2 | Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp) |

2291 | 1 | pfleura2 | Hess(IGeom,:,:)=HessTmp |

2292 | 1 | pfleura2 | END DO ! matches DO IGeom=2,NGeomF-1 |

2293 | 1 | pfleura2 | END IF !matches IF (HupNeighbour) THEN |

2294 | 1 | pfleura2 | END IF ! matches IF (FCalcHess) |

2295 | 1 | pfleura2 | END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc) |

2296 | 1 | pfleura2 | END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) |

2297 | 1 | pfleura2 | |

2298 | 2 | pfleura2 | ! IF (PROG=="GAUSSIAN") THEN |

2299 | 2 | pfleura2 | ! DEALLOCATE(Gauss_paste) |

2300 | 2 | pfleura2 | ! DEALLOCATE(Gauss_root) |

2301 | 2 | pfleura2 | ! DEALLOCATE(Gauss_comment) |

2302 | 2 | pfleura2 | ! DEALLOCATE(Gauss_end) |

2303 | 2 | pfleura2 | ! END IF |

2304 | 2 | pfleura2 | |

2305 | 1 | pfleura2 | DEALLOCATE(XyzGeomF, IntCoordF) |

2306 | 1 | pfleura2 | DEALLOCATE(XyzGeomI, IntCoordI) |

2307 | 1 | pfleura2 | DEALLOCATE(XyzTangent,IntTangent) |

2308 | 1 | pfleura2 | DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt) |

2309 | 1 | pfleura2 | DEALLOCATE(GeomOld_all) |

2310 | 1 | pfleura2 | |

2311 | 1 | pfleura2 | if (PROG=="GAUSSIAN") THEN |

2312 | 1 | pfleura2 | ! We de-allocate the Gauss_XX lists |

2313 | 1 | pfleura2 | ! transverse the list and deallocate each node |

2314 | 1 | pfleura2 | previous => Gauss_root ! make current point to head of list |

2315 | 1 | pfleura2 | DO |

2316 | 1 | pfleura2 | IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer |

2317 | 1 | pfleura2 | current => previous%next ! make list point to next node of head |

2318 | 1 | pfleura2 | DEALLOCATE(previous) ! deallocate current head node |

2319 | 1 | pfleura2 | previous => current ! make current point to new head |

2320 | 1 | pfleura2 | END DO |

2321 | 1 | pfleura2 | |

2322 | 1 | pfleura2 | previous => Gauss_comment ! make current point to head of list |

2323 | 1 | pfleura2 | DO |

2324 | 1 | pfleura2 | IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer |

2325 | 1 | pfleura2 | current => previous%next ! make list point to next node of head |

2326 | 1 | pfleura2 | DEALLOCATE(previous) ! deallocate current head node |

2327 | 1 | pfleura2 | previous => current ! make current point to new head |

2328 | 1 | pfleura2 | END DO |

2329 | 1 | pfleura2 | |

2330 | 1 | pfleura2 | previous => Gauss_end ! make current point to head of list |

2331 | 1 | pfleura2 | DO |

2332 | 1 | pfleura2 | IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer |

2333 | 1 | pfleura2 | current => previous%next ! make list point to next node of head |

2334 | 1 | pfleura2 | DEALLOCATE(previous) ! deallocate current head node |

2335 | 1 | pfleura2 | previous => current ! make current point to new head |

2336 | 1 | pfleura2 | END DO |

2337 | 1 | pfleura2 | |

2338 | 2 | pfleura2 | |

2339 | 1 | pfleura2 | END IF |

2340 | 1 | pfleura2 | |

2341 | 1 | pfleura2 | DEALLOCATE(GeomTmp, Hess, GradTmp) |

2342 | 1 | pfleura2 | IF (COORD.EQ.'ZMAT') DEALLOCATE(dzdc) |

2343 | 1 | pfleura2 | IF (COORD.EQ.'BAKER') THEN |

2344 | 1 | pfleura2 | DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t) |

2345 | 1 | pfleura2 | DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local) |

2346 | 1 | pfleura2 | END IF |

2347 | 1 | pfleura2 | |

2348 | 1 | pfleura2 | WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt |

2349 | 1 | pfleura2 | |

2350 | 7 | pfleura2 | Close(IOIN) |

2351 | 7 | pfleura2 | Close(IOOUT) |

2352 | 7 | pfleura2 | IF (AnaGeom) Close(IODAT) |

2353 | 8 | pfleura2 | IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot) |

2354 | 7 | pfleura2 | |

2355 | 1 | pfleura2 | END PROGRAM PathOpt |