Revision 8 src/Opt_Geom.f90
Opt_Geom.f90 (revision 8)  

3  3 
!so It has been designed to be almost independent of the rest of the code. 
4  4 
! It is now an (almost) officiel option. 
5  5  
6 
SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)


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

7  7  
8  8 
use VarTypes 
9  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 
use Io_module, only : IoGeom


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  14  
15 
use Io_module 

16  
15  17 
IMPLICIT NONE 
16  18  
17  19 
INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom 
...  ...  
30  32  
31  33 
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old) 
32  34  
33 
! This routines calculates the energy E and the gradient Grad of 

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

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

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.


36  38  
37 
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & 

38 
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, & 

39 
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & 

40 
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & 

41 
, Order,OrderInv, XPrimitiveF 

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


42  44  
43 
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) 

44 
! allocated in Path.f90 

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


46 
! allocated in Path.f90


45  47  
46 
use Io_module 

48 
use Io_module


47  49  
48 
! Energy (calculated if F300K=.F., else estimated) 

49 
REAL(KREAL), INTENT (OUT) :: E 

50 
! NCoord: Number of the degrees of freedom 

51 
! IGeom: index of the geometry. 

52 
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt 

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

54 
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) 

55 
! Gradient calculated at Geom geometry: 

56 
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) 

57 
! Cartesian geometry corresponding to (Internal Geometry) Geom: 

58 
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) 

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)


59  61 
!!! Optional, just for geometry optimization with Baker coordinates 
60 
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) 

61 
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) 

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

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

64 
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom 

65 
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F 

66 
LOGICAL :: Flag_Opt_Geom 

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


67  69  
68 
END subroutine Egrad 

70 
END subroutine Egrad


69  71  
70  72 
END INTERFACE 
71  73  
...  ...  
74  76 
LOGICAL :: Fini 
75  77 
LOGICAL, SAVE :: FRST=.TRUE. 
76  78  
77 
! Variables 

79 
! Variables


78  80 
INTEGER(KINT) :: IOpt, I,J,Iat, IBEG 
79  81 
REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat6 
80  82 
REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat6) 
...  ...  
89  91 
REAL(KREAL), ALLOCATABLE :: NewGradTmp(:) 
90  92 
REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis 
91  93 
REAL(KREAL) :: NormStep, FactStep, HP, MaxStep 
92 
REAL(KREAL) :: Eold 

94 
REAL(KREAL) :: Eold, Eini 

95 
! Values contains the values for the geometry analizes 

96 
REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar 

93  97  
94  98 
debug=valid('OptGeom') 
95  99  
96 
E=0. 

97 
Eold=0. 

98 
MaxStep=SMax 

100 
E=0. 

101 
Eold=0. 

102 
Eini=0. 

103 
MaxStep=SMax 

99  104  
100 
if (debug) Call Header("Entering Geom Optimization")


105 
if (debug) Call Header("Entering Geom Optimization") 

101  106  
102 
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord)) 

103 
ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord)) 

104 
ALLOCATE(GeomRef(NCoord)) 

105 
ALLOCATE(Hess_local(NCoord,NCoord)) 

106 
ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3)) 

107 
ALLOCATE(NewGeom(NCoord)) 

108 
ALLOCATE(NewGradTmp(NCoord)) 

109 
ALLOCATE(Hess_local_inv(NCoord,NCoord)) 

110 


111 
IOpt=0 

112 
ZeroVec=0.d0 

113 


114 
! Initialize the Hessian 

115 
Hess_local=0. 

107 
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord)) 

108 
ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord)) 

109 
ALLOCATE(GeomRef(NCoord)) 

110 
ALLOCATE(Hess_local(NCoord,NCoord)) 

111 
ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3)) 

112 
ALLOCATE(NewGeom(NCoord)) 

113 
ALLOCATE(NewGradTmp(NCoord)) 

114 
ALLOCATE(Hess_local_inv(NCoord,NCoord)) 

116  115  
117 
SELECT CASE (COORD) 

118 
CASE ('ZMAT') 

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

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

121 
Hess_local(1,1)=0.5d0 

122 
Hess_local(2,2)=0.5d0 

123 
Hess_local(3,3)=0.1d0 

124 
DO J=1,Nat3 

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

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

127 
Hess_local(3*J+3,3*J+3)=0.02d0 

116 
if (NbVar>0) THEN 

117 
ALLOCATE(Values(NbVar)) 

118 
END IF 

119 
FormAna(5:8)=" I5 " 

120 
IOpt=0 

121 
ZeroVec=0.d0 

122  
123 
! Initialize the Hessian 

124 
Hess_local=0. 

125  
126 
SELECT CASE (COORD) 

127 
CASE ('ZMAT') 

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

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

130 
Hess_local(1,1)=0.5d0 

131 
Hess_local(2,2)=0.5d0 

132 
Hess_local(3,3)=0.1d0 

133 
DO J=1,Nat3 

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

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

136 
Hess_local(3*J+3,3*J+3)=0.02d0 

137 
END DO 

138 
IF (HInv) THEN 

139 
DO I=1,NCoord 

140 
Hess_local(I,I)=1.d0/Hess_local(I,I) 

128  141 
END DO 
129 
IF (HInv) THEN 

130 
DO I=1,NCoord 

131 
Hess_local(I,I)=1.d0/Hess_local(I,I) 

132 
END DO 

133 
END IF 

142 
END IF 

134  143  
135 
IF (Step_method .EQ. "RFO") Then


136 
Hess_local=0.


137 
DO I=1,NCoord


138 
Hess_local(I,I)=0.5d0


139 
END DO


140 
END IF


144 
IF (Step_method .EQ. "RFO") Then 

145 
Hess_local=0. 

146 
DO I=1,NCoord 

147 
Hess_local(I,I)=0.5d0 

148 
END DO 

149 
END IF 

141  150  
142 
CASE ('BAKER') 

143 
! UMat(NPrim,3*Nat6) 

144 
BTransInv_local = BTransInv 

145 
UMat_local = UMat 

146 
ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord)) 

147 
Hprim=0.d0 

148 
ScanCoord=>Coordinate 

149 
I=0 

150 
DO WHILE (Associated(ScanCoord%next)) 

151 
I=I+1 

152 
SELECT CASE (ScanCoord%Type) 

153 
CASE ('BOND') 

154 
Hprim(I,I) = 0.5d0 

155 
CASE ('ANGLE') 

156 
Hprim(I,I) = 0.2d0 

157 
CASE ('DIHEDRAL') 

158 
Hprim(I,I) = 0.1d0 

159 
END SELECT 

160 
ScanCoord => ScanCoord%next 

161 
END DO 

162 


163 
! Hprim U: 

164 
HprimU=0.d0 

165 
DO I=1, NCoord 

166 
DO J=1,NPrim 

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

168 
END DO 

169 
END DO 

151 
CASE ('BAKER') 

152 
! UMat(NPrim,3*Nat6) 

153 
BTransInv_local = BTransInv 

154 
UMat_local = UMat 

155 
ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord)) 

156 
Hprim=0.d0 

157 
ScanCoord=>Coordinate 

158 
I=0 

159 
DO WHILE (Associated(ScanCoord%next)) 

160 
I=I+1 

161 
SELECT CASE (ScanCoord%Type) 

162 
CASE ('BOND') 

163 
Hprim(I,I) = 0.5d0 

164 
CASE ('ANGLE') 

165 
Hprim(I,I) = 0.2d0 

166 
CASE ('DIHEDRAL') 

167 
Hprim(I,I) = 0.1d0 

168 
END SELECT 

169 
ScanCoord => ScanCoord%next 

170 
END DO 

170  171  
171 
! Hess = U^T Hprim U: 

172 
Hess_local=0.d0 

173 
DO I=1, NCoord 

174 
DO J=1,NPrim 

175 
! UMat^T is needed below. 

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

177 
END DO 

178 
END DO 

179 


180 
!Print *, 'Hprim=' 

181 
DO I=1,NPrim 

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

183 
END DO 

184 
!Print *, 'UMat=' 

185 
DO I=1,NPrim 

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

187 
END DO 

188 
!Print *, 'HprimU=' 

189 
DO I=1,NPrim 

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

191 
END DO 

192 
!Print *, 'Hess_local=' 

193 
DO I=1,NCoord 

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

195 
END DO 

196 


197 
DEALLOCATE(Hprim,HprimU) 

198 


199 
IF (HInv) THEN 

200 
Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord) 

201 
Hess_local = Hess_local_inv 

202 
END IF 

203 


204 
!Print *, 'Hess_local after inversion=' 

205 
DO I=1,NCoord 

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

207 
END DO 

208 


209 
IF (Step_method .EQ. "RFO") Then 

210 
Hess_local=0. 

211 
DO I=1,NCoord 

212 
Hess_local(I,I)=0.5d0 

213 
END DO 

214 
END IF 

215 


216 
CASE ('MIXED','CART','HYBRID') 

217 
DO J=1,NCoord 

218 
Hess_local(J,J)=1. 

172 
! Hprim U: 

173 
HprimU=0.d0 

174 
DO I=1, NCoord 

175 
DO J=1,NPrim 

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

219  177 
END DO 
220 
CASE DEFAULT 

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

222 
STOP 

223 
END SELECT 

178 
END DO 

224  179  
225 
! Go to optimization 

226 
GeomOld=0.d0 ! Internal coordinates 

227 
GeomCart_old=0.d0 ! Cartesian coordinates 

180 
! Hess = U^T Hprim U: 

181 
Hess_local=0.d0 

182 
DO I=1, NCoord 

183 
DO J=1,NPrim 

184 
! UMat^T is needed below. 

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

186 
END DO 

187 
END DO 

228  188  
229 
IF (FPrintGeom) THEN 

230 
OPEN(IOGeom,File=GeomFile) 

189 
!Print *, 'Hprim=' 

190 
DO I=1,NPrim 

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

192 
END DO 

193 
!Print *, 'UMat=' 

194 
DO I=1,NPrim 

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

196 
END DO 

197 
!Print *, 'HprimU=' 

198 
DO I=1,NPrim 

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

200 
END DO 

201 
!Print *, 'Hess_local=' 

202 
DO I=1,NCoord 

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

204 
END DO 

205  
206 
DEALLOCATE(Hprim,HprimU) 

207  
208 
IF (HInv) THEN 

209 
Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord) 

210 
Hess_local = Hess_local_inv 

231  211 
END IF 
232  212  
233 
Fini=.FALSE. 

234 
DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) 

235 
IOpt=IOpt+1 

213 
!Print *, 'Hess_local after inversion=' 

214 
DO I=1,NCoord 

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

216 
END DO 

236  217  
237 
! Calculate the energy and gradient 

238 
IF (debug) THEN 

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

240 
WRITE(*,*) "Energy and Coordinates:" 

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

242 
WRITE(*,*) "Geom Old:" 

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

244 
WRITE(*,*) "Grad Old:" 

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

246 
END IF 

218 
IF (Step_method .EQ. "RFO") Then 

219 
Hess_local=0. 

220 
DO I=1,NCoord 

221 
Hess_local(I,I)=0.5d0 

222 
END DO 

223 
END IF 

247  224  
248 
!Call EGrad(E,Geom,GradTmp,NCoord) 

249 
IF (COORD.EQ.'BAKER') THEN 

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

251 
! GeomCart has INTENT(OUT) 

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

253 
GeomRef=GeomOld 

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

255 
GeomRef,GeomCart_old) 

256 
GeomCart_old=GeomCart 

225 
CASE ('MIXED','CART','HYBRID') 

226 
DO J=1,NCoord 

227 
Hess_local(J,J)=1. 

228 
END DO 

229 
CASE DEFAULT 

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

231 
STOP 

232 
END SELECT 

233  
234 
! Go to optimization 

235 
GeomOld=0.d0 ! Internal coordinates 

236 
GeomCart_old=0.d0 ! Cartesian coordinates 

237  
238 
IF (FPrintGeom) THEN 

239 
OPEN(IOGeom,File=GeomFile) 

240 
END IF 

241  
242 
Fini=.FALSE. 

243 
DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) 

244 
IOpt=IOpt+1 

245  
246 
! Calculate the energy and gradient 

247 
IF (debug) THEN 

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

249 
WRITE(*,*) "Energy and Coordinates:" 

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

251 
WRITE(*,*) "Geom Old:" 

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

253 
WRITE(*,*) "Grad Old:" 

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

255 
END IF 

256  
257 
!Call EGrad(E,Geom,GradTmp,NCoord) 

258 
IF (COORD.EQ.'BAKER') THEN 

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

260 
! GeomCart has INTENT(OUT) 

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

262 
GeomRef=GeomOld 

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

264 
GeomRef,GeomCart_old) 

265 
GeomCart_old=GeomCart 

266 
ELSE 

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

268 
END IF 

269  
270 
If (Iopt==1) EIni=E 

271  
272 
IF (debug) THEN 

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

274 
WRITE(*,*) "Energy and Coordinates:" 

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

276 
WRITE(*,*) "Geom Old:" 

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

278 
WRITE(*,*) "Grad:" 

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

280 
END IF 

281  
282 
IF (FPrintGeom) THEN 

283 
WRITE(IoGeom,*) Nat 

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

285  
286 
DO I=1,Nat 

287 
If (renum) THEN 

288 
Iat=Order(I) 

289 
WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat2:3*Iat) 

290 
ELSE 

291 
Iat=OrderInv(I) 

292 
WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I2:3*I) 

293 
END IF 

294 
END DO 

295 
END IF 

296  
297 
! do we have to analyze geometries ? 

298 
If (AnaGeom) THEN 

299 
If (NbVar>0) THEN 

300 
Call AnalyzeGeom(GeomCart,Values) 

301 
WRITE(IoDat,FormAna) Iopt,Values*PrintGeomFactor,(EEini)*ConvE,E 

302 
if (debug) WRITE(*,FormAna) Iopt,Values*PrintGeomFactor,(EEini)*ConvE,E 

257  303 
ELSE 
258 
Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart) 

304 
WRITE(IoDat,FormAna) Iopt,(EEini)*ConvE,E 

305 
if (debug) WRITE(*,FormAna) Iopt,(EEini)*ConvE,E 

259  306 
END IF 
260 


261 
IF (debug) THEN 

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

263 
WRITE(*,*) "Energy and Coordinates:" 

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

265 
WRITE(*,*) "Geom Old:" 

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

267 
WRITE(*,*) "Grad:" 

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

269 
END IF 

307 
END IF 

270  308  
271 
IF (FPrintGeom) THEN 

272 
WRITE(IoGeom,*) Nat 

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

274  309  
275 
DO I=1,Nat 

276 
If (renum) THEN 

277 
Iat=Order(I) 

278 
WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat2:3*Iat) 

279 
ELSE 

280 
Iat=OrderInv(I) 

281 
WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I2:3*I) 

282 
END IF 

283 
END DO 

284 
END IF 

285  
286 
IF (IOpt.GE.2) THEN 

310 
IF (IOpt.GE.2) THEN 

287  311 
! We have enough geometries and gradient to update the hessian (or its 
288  312 
! inverse) 
289  313 
GradTmp2=GradTmpGradOld 
...  ...  
291  315  
292  316 
if (debug) Write(*,*) "Call Hupdate_all, Geom" 
293  317 
IF (debug) THEN 
294 
WRITE(*,*) "dx before calling",SHAPE(GeomTmp2) 

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

296 
WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2) 

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

298 
END IF 

299 
Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local) 

300 
END IF ! matches IF (IOpt.GE.2) THEN 

301 


302 
GradOld=GradTmp 

303 
GeomOld=Geom 

318 
WRITE(*,*) "dx before calling",SHAPE(GeomTmp2) 

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

320 
WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2) 

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

322 
END IF 

323 
Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local) 

324 
END IF ! matches IF (IOpt.GE.2) THEN 

304  325  
305 
! GradTmp becomes Step in Step_RFO_all. 

306 
SELECT CASE (Step_method) 

307 
CASE ('RFO') 

308 
Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec) 

309 
CASE ('GDIIS') 

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

311 
! q_{m+1} = q'_{m+1}  H^{1}g'_{m+1} 

312 
Geom=0.d0 

313 
DO I=1, NCoord 

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

315 
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) 

316 
END DO 

317 
Geom(:) = NewGeom(:)  Geom(:) 

318 
! GradTmp now becomes "step" below: 

319 
GradTmp = Geom  GeomOld 

320 
CASE ('GDIIST') 

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

322 
! q_{m+1} = q'_{m+1}  H^{1}g'_{m+1} 

323 
Geom=0.d0 

324 
DO I=1, NCoord 

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

326 
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) 

327 
END DO 

328 
Geom(:) = NewGeom(:)  Geom(:) 

329 
! GradTmp now becomes "step" below: 

330 
GradTmp = Geom  GeomOld 

331 
CASE ('GEDIIS') 

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

333 
! GradTmp is actually "step" below: 

334 
GradTmp = NewGeom  GeomOld 

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

336 
CASE DEFAULT 

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

338 
STOP 

339 
END SELECT 

326 
GradOld=GradTmp 

327 
GeomOld=Geom 

340  328  
341 
Fini=.TRUE. 

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

343 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 

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

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

346 
Fini=(NormStep.LE.SThresh) 

347 
NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients. 

348 
WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh 

349 
Fini=(Fini.AND.(NormStep.LE.GThresh)) 

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

351 
GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all. 

329 
! GradTmp becomes Step in Step_RFO_all. 

330 
SELECT CASE (Step_method) 

331 
CASE ('RFO') 

332 
Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec) 

333 
CASE ('GDIIS') 

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

335 
! q_{m+1} = q'_{m+1}  H^{1}g'_{m+1} 

336 
Geom=0.d0 

337 
DO I=1, NCoord 

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

339 
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) 

340 
END DO 

341 
Geom(:) = NewGeom(:)  Geom(:) 

342 
! GradTmp now becomes "step" below: 

343 
GradTmp = Geom  GeomOld 

344 
CASE ('GDIIST') 

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

346 
! q_{m+1} = q'_{m+1}  H^{1}g'_{m+1} 

347 
Geom=0.d0 

348 
DO I=1, NCoord 

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

350 
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) 

351 
END DO 

352 
Geom(:) = NewGeom(:)  Geom(:) 

353 
! GradTmp now becomes "step" below: 

354 
GradTmp = Geom  GeomOld 

355 
CASE ('GEDIIS') 

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

357 
! GradTmp is actually "step" below: 

358 
GradTmp = NewGeom  GeomOld 

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

360 
CASE DEFAULT 

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

362 
STOP 

363 
END SELECT 

352  364  
365 
Fini=.TRUE. 

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

367 
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) 

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

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

370 
Fini=(NormStep.LE.SThresh) 

371 
NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients. 

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

373 
Fini=(Fini.AND.(NormStep.LE.GThresh)) 

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

375 
GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all. 

353  376  
354 
if (DynMaxStep.AND.(IOpt>1)) THEN 

355 
If (E<EOld) Then 

356 
MaxStep=min(1.1*MaxStep,2.*SMax) 

357 
ELSE 

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

359 
END IF 

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

361  377  
378 
if (DynMaxStep.AND.(IOpt>1)) THEN 

379 
If (E<EOld) Then 

380 
MaxStep=min(1.1*MaxStep,2.*SMax) 

381 
ELSE 

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

362  383 
END IF 
384 
if (debug) WRITE(*,*) "Dynamically updated Maximum step in Opt_Geom:",MaxStep,MaxStep/SMax 

363  385  
364 
Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp) 

365 


366 
Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all. 

386 
END IF 

367  387  
368 
EOld=E


388 
Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)


369  389  
370 
IF (debug) THEN 

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

372 
SELECT CASE (COORD) 

373 
CASE ('ZMAT','BAKER') 

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

375 
CASE('CART','HYBRID') 

376 
DO Iat=1,Nat 

377 
! PFL 30 Apr 2008 > 

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

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

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

381 
Geom(3*Iat2:3*Iat) 

382 
! < PFL 30 Apr 2008 

383 
END DO 

384 
CASE('MIXED') 

385 
DO Iat=1,NCart 

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

387 
Geom(3*Iat2:3*Iat) 

388 
END DO 

389 


390 
SELECT CASE (NCart) 

391 
CASE(1) 

392 
if (NAt.GE.2) THEN 

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

390 
Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all. 

391  
392 
EOld=E 

393  
394 
IF (debug) THEN 

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

396 
SELECT CASE (COORD) 

397 
CASE ('ZMAT','BAKER') 

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

399 
CASE('CART','HYBRID') 

400 
DO Iat=1,Nat 

401 
! PFL 30 Apr 2008 > 

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

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

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

405 
Geom(3*Iat2:3*Iat) 

406 
! < PFL 30 Apr 2008 

407 
END DO 

408 
CASE('MIXED') 

409 
DO Iat=1,NCart 

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

411 
Geom(3*Iat2:3*Iat) 

412 
END DO 

413  
414 
SELECT CASE (NCart) 

415 
CASE(1) 

416 
if (NAt.GE.2) THEN 

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

394  418 
IndZmat(2,2),Geom(4) 
395 
IBEG=5 

396 
END IF


397 
IF (NAT.GE.3) THEN


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

419 
IBEG=5


420 
END IF 

421 
IF (NAT.GE.3) THEN 

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


399  423 
IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6) 
400 
IBeg=7 

401 
END IF 

402 
CASE(2) 

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

404 
IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8) 

405 
IBeg=9 

406 
CASE DEFAULT 

407 
IBeg=3*NCart+1 

408 
END SELECT 

409  
410 
DO Iat=max(4,NCart),Nat 

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

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

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

414 
IBeg=IBeg+3 

415 
END DO 

416  
424 
IBeg=7 

425 
END IF 

426 
CASE(2) 

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

428 
IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8) 

429 
IBeg=9 

417  430 
CASE DEFAULT 
418 
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294." 

419 
STOP 

431 
IBeg=3*NCart+1 

420  432 
END SELECT 
421 
END IF ! matches IF (debug) THEN 

422 


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

424 


425 
DEALLOCATE(GradTmp2,GeomTmp2,GradTmp) 

426 
DEALLOCATE(GradOld,GeomOld) 

427 
DEALLOCATE(Hess_local) 

428 
DEALLOCATE(GeomCart_old) 

429 
DEALLOCATE(NewGeom,NewGradTmp) 

430 
DEALLOCATE(Hess_local_inv) 

431  433  
432 
IF (Fini) THEN 

433 
WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E 

434 
ELSE 

435 
WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E 

436 
END IF 

437 


438 
WRITE(*,*) "Initial Geometry:" 

439 
DO I=1, Nat 

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

441 
END DO 

442 
WRITE(*,*) "Final Geometry:" 

443 
DO I=1, Nat 

444 
WRITE(*,'(3(1X,F8.4))') GeomCart(I,:) 

445 
!IF (I .GT. 1) Then 

446 
! WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)GeomCart(I1,1))*(GeomCart(I,1)GeomCart(I1,1)))& 

447 
! + ((GeomCart(I,2)GeomCart(I1,2))*(GeomCart(I,2)GeomCart(I1,2)))& 

448 
! + ((GeomCart(I,3)GeomCart(I1,3))*(GeomCart(I,3)GeomCart(I1,3)))) 

449 
!END IF 

450 
END DO 

451 


452 
IF (COORD .EQ. "BAKER") Then 

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

454 
ScanCoord=>Coordinate 

455 
DO WHILE (Associated(ScanCoord%next)) 

456 
I=I+1 

457 
SELECT CASE (ScanCoord%Type) 

458 
CASE ('BOND') 

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

460 
CASE ('ANGLE') 

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

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

463 
CASE ('DIHEDRAL') 

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

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

466 
END SELECT 

467 
ScanCoord => ScanCoord%next 

468 
END DO ! matches DO WHILE (Associated(ScanCoord%next)) 

469 
END IF 

470 


471 
DEALLOCATE(GeomCart) 

434 
DO Iat=max(4,NCart),Nat 

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

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

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

438 
IBeg=IBeg+3 

439 
END DO 

472  440  
473 
if (debug) Call Header("Geom Optimization Over") 

441 
CASE DEFAULT 

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

443 
STOP 

444 
END SELECT 

445 
END IF ! matches IF (debug) THEN 

474  446  
475 
END SUBROUTINE Opt_geom 

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

448  
449 
DEALLOCATE(GradTmp2,GeomTmp2,GradTmp) 

450 
DEALLOCATE(GradOld,GeomOld) 

451 
DEALLOCATE(Hess_local) 

452 
DEALLOCATE(GeomCart_old) 

453 
DEALLOCATE(NewGeom,NewGradTmp) 

454 
DEALLOCATE(Hess_local_inv) 

455  
456 
IF (Fini) THEN 

457 
WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E 

458 
ELSE 

459 
WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E 

460 
END IF 

461  
462 
WRITE(*,*) "Initial Geometry:" 

463 
DO I=1, Nat 

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

465 
END DO 

466 
WRITE(*,*) "Final Geometry:" 

467 
DO I=1, Nat 

468 
WRITE(*,'(3(1X,F8.4))') GeomCart(I,:) 

469 
!IF (I .GT. 1) Then 

470 
! WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)GeomCart(I1,1))*(GeomCart(I,1)GeomCart(I1,1)))& 

471 
! + ((GeomCart(I,2)GeomCart(I1,2))*(GeomCart(I,2)GeomCart(I1,2)))& 

472 
! + ((GeomCart(I,3)GeomCart(I1,3))*(GeomCart(I,3)GeomCart(I1,3)))) 

473 
!END IF 

474 
END DO 

475  
476 
IF (COORD .EQ. "BAKER") Then 

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

478 
ScanCoord=>Coordinate 

479 
DO WHILE (Associated(ScanCoord%next)) 

480 
I=I+1 

481 
SELECT CASE (ScanCoord%Type) 

482 
CASE ('BOND') 

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

484 
CASE ('ANGLE') 

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

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

487 
CASE ('DIHEDRAL') 

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

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

490 
END SELECT 

491 
ScanCoord => ScanCoord%next 

492 
END DO ! matches DO WHILE (Associated(ScanCoord%next)) 

493 
END IF 

494  
495 
DEALLOCATE(GeomCart) 

496  
497 
if (debug) Call Header("Geom Optimization Over") 

498  
499 
END SUBROUTINE Opt_geom 
Also available in: Unified diff