## root / src / Opt_Geom.f90 @ 9

Historique | Voir | Annoter | Télécharger (18,94 ko)

1 |
! This subroutine optimizes a geometry. |
---|---|

2 |
! Initially, it was mainly here for debugging purposes.. |

3 |
!so It has been designed to be almost independent of the rest of the code. |

4 |
! It is now an (almost) officiel option. |

5 | |

6 |
SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method) |

7 | |

8 |
use VarTypes |

9 |
use Path_module, only : maxcyc,sthresh,gthresh,Hinv,Smax,Nom,Atome,OrderInv,indzmat,& |

10 |
XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,& |

11 |
BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,& |

12 |
Coordinate,Pi,UMat_local,NCart, DynMaxStep,FPrintGeom,GeomFile, AtName,Renum,Order, & |

13 |
FormAna,NbVar, PrintGeomFactor,AnaGeom |

14 | |

15 |
use Io_module |

16 | |

17 |
IMPLICIT NONE |

18 | |

19 |
INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom |

20 |
CHARACTER(32), INTENT(IN) :: Coord |

21 |
CHARACTER(32), INTENT(IN) :: Step_method |

22 |
REAL(KREAL), INTENT(INOUT) :: Geom(NCoord) |

23 |
REAL(KREAL), INTENT(OUT) :: E |

24 |
LOGICAL, INTENT(IN) :: Flag_Opt_Geom |

25 | |

26 | |

27 |
INTERFACE |

28 |
function valid(string) result (isValid) |

29 |
CHARACTER(*), intent(in) :: string |

30 |
logical :: isValid |

31 |
END function VALID |

32 | |

33 |
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old) |

34 | |

35 |
! This routines calculates the energy E and the gradient Grad of |

36 |
! a molecule with Geometry Geom (may be in internal coordinates), |

37 |
! using for now, either Gaussian or Ext, more general later. |

38 | |

39 |
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & |

40 |
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, & |

41 |
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & |

42 |
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & |

43 |
, Order,OrderInv, XPrimitiveF |

44 | |

45 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |

46 |
! allocated in Path.f90 |

47 | |

48 |
use Io_module |

49 | |

50 |
! Energy (calculated if F300K=.F., else estimated) |

51 |
REAL(KREAL), INTENT (OUT) :: E |

52 |
! NCoord: Number of the degrees of freedom |

53 |
! IGeom: index of the geometry. |

54 |
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt |

55 |
! Geometry at which gradient is calculated (cf Factual also): |

56 |
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) |

57 |
! Gradient calculated at Geom geometry: |

58 |
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) |

59 |
! Cartesian geometry corresponding to (Internal Geometry) Geom: |

60 |
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) |

61 |
!!! Optional, just for geometry optimization with Baker coordinates |

62 |
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) |

63 |
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) |

64 |
! FOptGeom is a flag indicating if we are doing a geom optimization |

65 |
! it can be omitted so that we use a local flag: Flag_Opt_Geom |

66 |
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom |

67 |
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F |

68 |
LOGICAL :: Flag_Opt_Geom |

69 | |

70 |
END subroutine Egrad |

71 | |

72 |
END INTERFACE |

73 | |

74 | |

75 |
LOGICAL :: debug |

76 |
LOGICAL :: FiniS,FiniG,Fini |

77 |
LOGICAL, SAVE :: FRST=.TRUE. |

78 | |

79 |
! Variables |

80 |
INTEGER(KINT) :: IOpt, I,J,Iat, IBEG |

81 |
REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6 |

82 |
REAL(KREAL), ALLOCATABLE :: Step(:) ! 3*Nat or 3*Nat-6 |

83 |
REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6) |

84 |
REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6) |

85 |
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6 |

86 |
REAL(KREAL), ALLOCATABLE :: Hess_local(:,:) |

87 |
REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3) |

88 |
REAL(KREAL), ALLOCATABLE :: GeomCart_old(:,:) !(Nat,3) |

89 |
REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim) |

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

91 |
REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6 |

92 |
REAL(KREAL), ALLOCATABLE :: NewGradTmp(:) |

93 |
REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis |

94 |
REAL(KREAL) :: NormStep, FactStep, HP, MaxStep,NormGrad |

95 |
REAL(KREAL) :: Eold, Eini |

96 |
! Values contains the values for the geometry analizes |

97 |
REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar |

98 |
CHARACTER(LCHARS) :: Line |

99 | |

100 |
debug=valid('OptGeom') |

101 | |

102 |
E=0. |

103 |
Eold=0. |

104 |
Eini=0. |

105 |
MaxStep=SMax |

106 | |

107 |
if (debug) Call Header("Entering Geom Optimization") |

108 | |

109 |
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord),Step(NCoord)) |

110 |
ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord)) |

111 |
ALLOCATE(GeomRef(NCoord)) |

112 |
ALLOCATE(Hess_local(NCoord,NCoord)) |

113 |
ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3)) |

114 |
ALLOCATE(NewGeom(NCoord)) |

115 |
ALLOCATE(NewGradTmp(NCoord)) |

116 |
ALLOCATE(Hess_local_inv(NCoord,NCoord)) |

117 | |

118 |
if (NbVar>0) THEN |

119 |
ALLOCATE(Values(NbVar)) |

120 |
END IF |

121 |
FormAna(5:8)=" I5 " |

122 |
IOpt=0 |

123 |
ZeroVec=0.d0 |

124 | |

125 |
! Initialize the Hessian |

126 |
Hess_local=0. |

127 | |

128 |
SELECT CASE (COORD) |

129 |
CASE ('ZMAT') |

130 |
! We use the fact that we know that approximate good hessian values |

131 |
! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles |

132 |
Hess_local(1,1)=0.5d0 |

133 |
Hess_local(2,2)=0.5d0 |

134 |
Hess_local(3,3)=0.1d0 |

135 |
DO J=1,Nat-3 |

136 |
Hess_local(3*J+1,3*J+1)=0.5d0 |

137 |
Hess_local(3*J+2,3*J+2)=0.1d0 |

138 |
Hess_local(3*J+3,3*J+3)=0.02d0 |

139 |
END DO |

140 |
IF (HInv) THEN |

141 |
DO I=1,NCoord |

142 |
Hess_local(I,I)=1.d0/Hess_local(I,I) |

143 |
END DO |

144 |
END IF |

145 | |

146 |
IF (Step_method .EQ. "RFO") Then |

147 |
Hess_local=0. |

148 |
DO I=1,NCoord |

149 |
Hess_local(I,I)=1.d0 |

150 |
END DO |

151 |
END IF |

152 | |

153 |
CASE ('BAKER') |

154 |
! UMat(NPrim,3*Nat-6) |

155 |
BTransInv_local = BTransInv |

156 |
UMat_local = UMat |

157 |
ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord)) |

158 |
Hprim=0.d0 |

159 |
ScanCoord=>Coordinate |

160 |
I=0 |

161 |
DO WHILE (Associated(ScanCoord%next)) |

162 |
I=I+1 |

163 |
SELECT CASE (ScanCoord%Type) |

164 |
CASE ('BOND') |

165 |
Hprim(I,I) = 0.5d0 |

166 |
CASE ('ANGLE') |

167 |
Hprim(I,I) = 0.2d0 |

168 |
CASE ('DIHEDRAL') |

169 |
Hprim(I,I) = 0.1d0 |

170 |
END SELECT |

171 |
ScanCoord => ScanCoord%next |

172 |
END DO |

173 | |

174 |
! Hprim U: |

175 |
HprimU=0.d0 |

176 |
DO I=1, NCoord |

177 |
DO J=1,NPrim |

178 |
HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I) |

179 |
END DO |

180 |
END DO |

181 | |

182 |
! Hess = U^T Hprim U: |

183 |
Hess_local=0.d0 |

184 |
DO I=1, NCoord |

185 |
DO J=1,NPrim |

186 |
! UMat^T is needed below. |

187 |
Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I) |

188 |
END DO |

189 |
END DO |

190 | |

191 |
!Print *, 'Hprim=' |

192 |
DO I=1,NPrim |

193 |
! WRITE(*,'(10(1X,F6.3))') Hprim(I,:) |

194 |
END DO |

195 |
!Print *, 'UMat=' |

196 |
DO I=1,NPrim |

197 |
! WRITE(*,'(3(1X,F6.3))') UMat(I,1:3) |

198 |
END DO |

199 |
!Print *, 'HprimU=' |

200 |
DO I=1,NPrim |

201 |
! WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3) |

202 |
END DO |

203 |
!Print *, 'Hess_local=' |

204 |
DO I=1,NCoord |

205 |
! WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3) |

206 |
END DO |

207 | |

208 |
DEALLOCATE(Hprim,HprimU) |

209 | |

210 |
IF (HInv) THEN |

211 |
Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord) |

212 |
Hess_local = Hess_local_inv |

213 |
END IF |

214 | |

215 |
!Print *, 'Hess_local after inversion=' |

216 |
! DO I=1,NCoord |

217 |
! WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3) |

218 |
! END DO |

219 | |

220 |
IF (Step_method .EQ. "RFO") Then |

221 |
Hess_local=0. |

222 |
DO I=1,NCoord |

223 |
Hess_local(I,I)=0.5d0 |

224 |
END DO |

225 |
END IF |

226 | |

227 |
CASE ('MIXED','CART','HYBRID') |

228 |
DO J=1,NCoord |

229 |
Hess_local(J,J)=1. |

230 |
END DO |

231 |
CASE DEFAULT |

232 |
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop" |

233 |
STOP |

234 |
END SELECT |

235 | |

236 |
! Go to optimization |

237 |
GeomOld=0.d0 ! Internal coordinates |

238 |
GeomCart_old=0.d0 ! Cartesian coordinates |

239 | |

240 |
IF (FPrintGeom) THEN |

241 |
OPEN(IOGeom,File=GeomFile) |

242 |
END IF |

243 | |

244 |
Fini=.FALSE. |

245 |
DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) |

246 |
IOpt=IOpt+1 |

247 | |

248 |
Write(Line,'(1X,A,I5)') "Iteration ",Iopt |

249 |
Call Header(TRIM(Line)) |

250 |
WRITE(IoOut,*) "Current Geometry" |

251 |
Line="" |

252 |
Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat) |

253 |
if (COORD/="CART") THEN |

254 |
WRITE(IoOut,*) "Current Geometry in Cartesian coordinates" |

255 |
Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat) |

256 |
END IF |

257 | |

258 |
WRITE(IoOut,*) "Computing energy and gradient" |

259 |
! Calculate the energy and gradient |

260 |
IF (debug) THEN |

261 |
WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt |

262 |
WRITE(*,*) "Energy and Coordinates:" |

263 |
WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12)) |

264 |
WRITE(*,*) "Geom Old:" |

265 |
WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13)) |

266 |
WRITE(*,*) "Grad Old:" |

267 |
WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13)) |

268 |
END IF |

269 | |

270 |
!Call EGrad(E,Geom,GradTmp,NCoord) |

271 |
IF (COORD.EQ.'BAKER') THEN |

272 |
! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT). |

273 |
! GeomCart has INTENT(OUT) |

274 |
! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable |

275 |
GeomRef=GeomOld |

276 |
Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, & |

277 |
GeomRef,GeomCart_old) |

278 |
GeomCart_old=GeomCart |

279 |
ELSE |

280 |
Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart) |

281 |
END IF |

282 |
!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!! |

283 |
! We have a pb with the format of GradTmp(Nat,3) whereas it is 3,Nat in Egrad |

284 |
! |

285 |
! This is a path for CART coordinates !!! |

286 |
IF (COORD=="CART") THEN |

287 |
Step=reshape(reshape(GradTmp,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/)) |

288 |
GradTmp=Step |

289 |
END IF |

290 |
!!!!!!!!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!!!!!!!! |

291 | |

292 |
WRITE(IoOut,'(1X," Energy E=",F15.6,A)') E,TRIM(UnitProg) |

293 |
WRITE(IoOut,*) "Gradient for current geometry" |

294 |
Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,IOOUT,Atome,Order,OrderInv,IndZmat) |

295 | |

296 |
If (Iopt==1) EIni=E |

297 | |

298 |
IF (debug) THEN |

299 |
WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt |

300 |
WRITE(*,*) "Energy and Coordinates:" |

301 |
WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12)) |

302 |
WRITE(*,*) "Geom Old:" |

303 |
WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12)) |

304 |
WRITE(*,*) "Grad:" |

305 |
WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12)) |

306 |
END IF |

307 | |

308 |
IF (FPrintGeom) THEN |

309 |
WRITE(IoGeom,*) Nat |

310 |
WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E |

311 | |

312 |
DO I=1,Nat |

313 |
If (renum) THEN |

314 |
Iat=Order(I) |

315 |
WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat) |

316 |
ELSE |

317 |
Iat=OrderInv(I) |

318 |
WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I) |

319 |
END IF |

320 |
END DO |

321 |
END IF |

322 | |

323 |
! do we have to analyze geometries ? |

324 |
If (AnaGeom) THEN |

325 |
If (NbVar>0) THEN |

326 |
Call AnalyzeGeom(GeomCart,Values) |

327 |
WRITE(IoDat,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E |

328 |
if (debug) WRITE(*,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E |

329 |
ELSE |

330 |
WRITE(IoDat,FormAna) Iopt,(E-Eini)*ConvE,E |

331 |
if (debug) WRITE(*,FormAna) Iopt,(E-Eini)*ConvE,E |

332 |
END IF |

333 |
END IF |

334 | |

335 | |

336 |
IF (IOpt.GE.2) THEN |

337 |
! We have enough geometries and gradient to update the hessian (or its |

338 |
! inverse) |

339 |
GradTmp2=GradTmp-GradOld |

340 |
GeomTmp2=Geom-GeomOld |

341 | |

342 |
if (debug) Write(*,*) "Call Hupdate_all, Geom" |

343 |
IF (debug) THEN |

344 |
WRITE(*,*) "dx before calling",SHAPE(GeomTmp2) |

345 |
WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord) |

346 |
WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2) |

347 |
WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord) |

348 |
END IF |

349 |
Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local) |

350 |
END IF ! matches IF (IOpt.GE.2) THEN |

351 | |

352 |
GradOld=GradTmp |

353 |
GeomOld=Geom |

354 | |

355 |
! GradTmp becomes Step in Step_RFO_all. |

356 |
SELECT CASE (Step_method) |

357 |
CASE ('RFO') |

358 |
Call Step_RFO_all(NCoord,Step,IGeom,Geom,GradTmp,Hess_local,ZeroVec) |

359 |
GradTmp=Step |

360 |
CASE ('GDIIS') |

361 |
Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST) |

362 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1} |

363 |
Geom=0.d0 |

364 |
DO I=1, NCoord |

365 |
! If Hinv=.False., then we need to invert Hess_local |

366 |
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) |

367 |
END DO |

368 |
Geom(:) = NewGeom(:) - Geom(:) |

369 |
! GradTmp now becomes "step" below: |

370 |
GradTmp = Geom - GeomOld |

371 |
CASE ('GDIIST') |

372 |
Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,NCoord,FRST) |

373 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1} |

374 |
Geom=0.d0 |

375 |
DO I=1, NCoord |

376 |
! If Hinv=.False., then we need to invert Hess_local |

377 |
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) |

378 |
END DO |

379 |
Geom(:) = NewGeom(:) - Geom(:) |

380 |
! GradTmp now becomes "step" below: |

381 |
GradTmp = Geom - GeomOld |

382 |
CASE ('GEDIIS') |

383 |
Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST) |

384 |
! GradTmp is actually "step" below: |

385 |
GradTmp = NewGeom - GeomOld |

386 |
!Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec) |

387 |
CASE DEFAULT |

388 |
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop" |

389 |
STOP |

390 |
END SELECT |

391 | |

392 |
Fini=.TRUE. |

393 |
! If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0 |

394 |
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) |

395 |
if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh |

396 |
FactStep=min(1.d0,MaxStep/NormStep) |

397 |
FiniS=((NormStep*FactStep).LE.SThresh) |

398 |
NormGrad=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients. |

399 |
if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormGrad (Gradient)=', NormGrad, 'Gradient Threshold=', GThresh |

400 |
FiniG=(NormGrad.LE.GThresh) |

401 |
Fini=FiniS.AND.FiniG |

402 |
if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep |

403 |
GradTmp=GradTmp*FactStep ! GradTmp is step here, from Step_RFO_all. |

404 | |

405 |
WRITE(IoOut,*) " Converged ?" |

406 |
WRITE(IoOut,'(1X,A18,1X,A15,1X,A13,5X,A3)') " Criterion","Current Value","Threshold","Ok?" |

407 |
IF (FiniS) THEN |

408 |
WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"YES" |

409 |
ELSE |

410 |
WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"NO" |

411 |
END IF |

412 |
IF (FiniG) THEN |

413 |
WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"YES" |

414 |
ELSE |

415 |
WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"NO" |

416 |
END IF |

417 | |

418 |
if (DynMaxStep.AND.(IOpt>1)) THEN |

419 |
If (E<EOld) Then |

420 |
MaxStep=min(1.1*MaxStep,2.*SMax) |

421 |
ELSE |

422 |
MaxStep=max(0.8*MaxStep,SMax/2.) |

423 |
END IF |

424 |
if (debug) WRITE(*,*) "Dynamically updated Maximum step in Opt_Geom:",MaxStep,MaxStep/SMax |

425 | |

426 |
END IF |

427 | |

428 |
Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp) |

429 | |

430 |
! We add the step to the old geom |

431 |
Geom=GeomOld+GradTmp |

432 | |

433 |
EOld=E |

434 | |

435 |
IF (debug) THEN |

436 |
WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt |

437 |
SELECT CASE (COORD) |

438 |
CASE ('ZMAT','BAKER') |

439 |
!WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord)) |

440 |
CASE('CART','HYBRID') |

441 |
DO Iat=1,Nat |

442 |
! PFL 30 Apr 2008 -> |

443 |
! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid. |

444 |
! WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |

445 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), & |

446 |
Geom(3*Iat-2:3*Iat) |

447 |
! <- PFL 30 Apr 2008 |

448 |
END DO |

449 |
CASE('MIXED') |

450 |
DO Iat=1,NCart |

451 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |

452 |
Geom(3*Iat-2:3*Iat) |

453 |
END DO |

454 | |

455 |
SELECT CASE (NCart) |

456 |
CASE(1) |

457 |
if (NAt.GE.2) THEN |

458 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & |

459 |
IndZmat(2,2),Geom(4) |

460 |
IBEG=5 |

461 |
END IF |

462 |
IF (NAT.GE.3) THEN |

463 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |

464 |
IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6) |

465 |
IBeg=7 |

466 |
END IF |

467 |
CASE(2) |

468 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |

469 |
IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8) |

470 |
IBeg=9 |

471 |
CASE DEFAULT |

472 |
IBeg=3*NCart+1 |

473 |
END SELECT |

474 | |

475 |
DO Iat=max(4,NCart),Nat |

476 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |

477 |
IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), & |

478 |
IndZmat(3,3),Geom(IBeg+2)*180./pi |

479 |
IBeg=IBeg+3 |

480 |
END DO |

481 | |

482 |
CASE DEFAULT |

483 |
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294." |

484 |
STOP |

485 |
END SELECT |

486 |
END IF ! matches IF (debug) THEN |

487 | |

488 |
END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) |

489 | |

490 |
IF (Fini) THEN |

491 |
WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles" |

492 |
ELSE |

493 |
WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," *NOT* converged in ",Iopt," cycles" |

494 |
END IF |

495 |
Call Header(Line) |

496 |
WRITE(IoOut,*) "Last Energy E=",E |

497 | |

498 |
! WRITE(*,*) "Initial Geometry:" |

499 |
! DO I=1, Nat |

500 |
! WRITE(*,'(3(1X,F8.4))') XyzGeomI(IGeom,:,I) |

501 |
! END DO |

502 |
WRITE(IoOut,*) "Final Geometry:" |

503 | |

504 |
Line="" |

505 |
Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat) |

506 |
if (COORD/="CART") THEN |

507 |
WRITE(IoOut,*) "Last Geometry in Cartesian coordinates" |

508 |
Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat) |

509 |
END IF |

510 | |

511 | |

512 |
! DO I=1, Nat |

513 |
! WRITE(*,'(3(1X,F8.4))') GeomCart(I,:) |

514 |
!IF (I .GT. 1) Then |

515 |
! WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))& |

516 |
! + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))& |

517 |
! + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3)))) |

518 |
!END IF |

519 |
! END DO |

520 | |

521 |
IF (COORD .EQ. "BAKER") Then |

522 |
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |

523 |
ScanCoord=>Coordinate |

524 |
DO WHILE (Associated(ScanCoord%next)) |

525 |
I=I+1 |

526 |
SELECT CASE (ScanCoord%Type) |

527 |
CASE ('BOND') |

528 |
WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I) |

529 |
CASE ('ANGLE') |

530 |
WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, & |

531 |
ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi |

532 |
CASE ('DIHEDRAL') |

533 |
WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,& |

534 |
ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi |

535 |
END SELECT |

536 |
ScanCoord => ScanCoord%next |

537 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |

538 |
END IF |

539 | |

540 |
DEALLOCATE(GeomCart) |

541 |
DEALLOCATE(GradTmp2,GeomTmp2,GradTmp,Step) |

542 |
DEALLOCATE(GradOld,GeomOld) |

543 |
DEALLOCATE(Hess_local) |

544 |
DEALLOCATE(GeomCart_old) |

545 |
DEALLOCATE(NewGeom,NewGradTmp) |

546 |
DEALLOCATE(Hess_local_inv) |

547 | |

548 |
if (debug) Call Header("Geom Optimization Over") |

549 | |

550 |
END SUBROUTINE Opt_geom |