root / src / Opt_Geom.f90 @ 5
Historique | Voir | Annoter | Télécharger (14,95 ko)
1 |
! This subroutine optimizes a geometry. It is mainly here for debugging purposes.. |
---|---|
2 |
! It has been designed to be almost independent of the rest of the code. |
3 |
|
4 |
SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method) |
5 |
|
6 |
use Path_module, only : maxcyc,sthresh,gthresh,Hinv,Smax,Nom,Atome,OrderInv,indzmat,& |
7 |
XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,& |
8 |
BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,& |
9 |
Coordinate,Pi,UMat_local,NCart |
10 |
use VarTypes |
11 |
IMPLICIT NONE |
12 |
|
13 |
INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom |
14 |
CHARACTER(32), INTENT(IN) :: Coord |
15 |
CHARACTER(32), INTENT(IN) :: Step_method |
16 |
REAL(KREAL), INTENT(INOUT) :: Geom(NCoord) |
17 |
REAL(KREAL), INTENT(OUT) :: E |
18 |
LOGICAL, INTENT(IN) :: Flag_Opt_Geom |
19 |
|
20 |
! As this subroutine is here for debugging, debug=T ! |
21 |
LOGICAL, PARAMETER :: debug=.TRUE. |
22 |
LOGICAL :: Fini |
23 |
LOGICAL, SAVE :: FRST=.TRUE. |
24 |
|
25 |
! Variables |
26 |
INTEGER(KINT) :: IOpt, I,J,Iat, IBEG |
27 |
REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6 |
28 |
REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6) |
29 |
REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6) |
30 |
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6 |
31 |
REAL(KREAL), ALLOCATABLE :: Hess_local(:,:) |
32 |
REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3) |
33 |
REAL(KREAL), ALLOCATABLE :: GeomCart_old(:,:) !(Nat,3) |
34 |
REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim) |
35 |
REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6)) |
36 |
REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6 |
37 |
REAL(KREAL), ALLOCATABLE :: NewGradTmp(:) |
38 |
REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis |
39 |
REAL(KREAL) :: NormStep, FactStep, HP |
40 |
|
41 |
E=0. |
42 |
|
43 |
if (debug) WRITE(*,*) "=================== Geom Optimization ==================" |
44 |
|
45 |
ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord)) |
46 |
ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord)) |
47 |
ALLOCATE(GeomRef(NCoord)) |
48 |
ALLOCATE(Hess_local(NCoord,NCoord)) |
49 |
ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3)) |
50 |
ALLOCATE(NewGeom(NCoord)) |
51 |
ALLOCATE(NewGradTmp(NCoord)) |
52 |
ALLOCATE(Hess_local_inv(NCoord,NCoord)) |
53 |
|
54 |
IOpt=0 |
55 |
ZeroVec=0.d0 |
56 |
|
57 |
! Initialize the Hessian |
58 |
Hess_local=0. |
59 |
|
60 |
SELECT CASE (COORD) |
61 |
CASE ('ZMAT') |
62 |
! We use the fact that we know that approximate good hessian values |
63 |
! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles |
64 |
Hess_local(1,1)=0.5d0 |
65 |
Hess_local(2,2)=0.5d0 |
66 |
Hess_local(3,3)=0.1d0 |
67 |
DO J=1,Nat-3 |
68 |
Hess_local(3*J+1,3*J+1)=0.5d0 |
69 |
Hess_local(3*J+2,3*J+2)=0.1d0 |
70 |
Hess_local(3*J+3,3*J+3)=0.02d0 |
71 |
END DO |
72 |
IF (HInv) THEN |
73 |
DO I=1,NCoord |
74 |
Hess_local(I,I)=1.d0/Hess_local(I,I) |
75 |
END DO |
76 |
END IF |
77 |
|
78 |
IF (Step_method .EQ. "RFO") Then |
79 |
Hess_local=0. |
80 |
DO I=1,NCoord |
81 |
Hess_local(I,I)=0.5d0 |
82 |
END DO |
83 |
END IF |
84 |
|
85 |
CASE ('BAKER') |
86 |
! UMat(NPrim,3*Nat-6) |
87 |
BTransInv_local = BTransInv |
88 |
UMat_local = UMat |
89 |
ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord)) |
90 |
Hprim=0.d0 |
91 |
ScanCoord=>Coordinate |
92 |
I=0 |
93 |
DO WHILE (Associated(ScanCoord%next)) |
94 |
I=I+1 |
95 |
SELECT CASE (ScanCoord%Type) |
96 |
CASE ('BOND') |
97 |
Hprim(I,I) = 0.5d0 |
98 |
CASE ('ANGLE') |
99 |
Hprim(I,I) = 0.2d0 |
100 |
CASE ('DIHEDRAL') |
101 |
Hprim(I,I) = 0.1d0 |
102 |
END SELECT |
103 |
ScanCoord => ScanCoord%next |
104 |
END DO |
105 |
|
106 |
! Hprim U: |
107 |
HprimU=0.d0 |
108 |
DO I=1, NCoord |
109 |
DO J=1,NPrim |
110 |
HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I) |
111 |
END DO |
112 |
END DO |
113 |
|
114 |
! Hess = U^T Hprim U: |
115 |
Hess_local=0.d0 |
116 |
DO I=1, NCoord |
117 |
DO J=1,NPrim |
118 |
! UMat^T is needed below. |
119 |
Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I) |
120 |
END DO |
121 |
END DO |
122 |
|
123 |
!Print *, 'Hprim=' |
124 |
DO I=1,NPrim |
125 |
! WRITE(*,'(10(1X,F6.3))') Hprim(I,:) |
126 |
END DO |
127 |
!Print *, 'UMat=' |
128 |
DO I=1,NPrim |
129 |
! WRITE(*,'(3(1X,F6.3))') UMat(I,1:3) |
130 |
END DO |
131 |
!Print *, 'HprimU=' |
132 |
DO I=1,NPrim |
133 |
! WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3) |
134 |
END DO |
135 |
!Print *, 'Hess_local=' |
136 |
DO I=1,NCoord |
137 |
! WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3) |
138 |
END DO |
139 |
|
140 |
DEALLOCATE(Hprim,HprimU) |
141 |
|
142 |
IF (HInv) THEN |
143 |
Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord) |
144 |
Hess_local = Hess_local_inv |
145 |
END IF |
146 |
|
147 |
!Print *, 'Hess_local after inversion=' |
148 |
DO I=1,NCoord |
149 |
! WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3) |
150 |
END DO |
151 |
|
152 |
IF (Step_method .EQ. "RFO") Then |
153 |
Hess_local=0. |
154 |
DO I=1,NCoord |
155 |
Hess_local(I,I)=0.5d0 |
156 |
END DO |
157 |
END IF |
158 |
|
159 |
CASE ('MIXED','CART','HYBRID') |
160 |
DO J=1,NCoord |
161 |
Hess_local(J,J)=1. |
162 |
END DO |
163 |
CASE DEFAULT |
164 |
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop" |
165 |
STOP |
166 |
END SELECT |
167 |
|
168 |
! Go to optimization |
169 |
GeomOld=0.d0 ! Internal coordinates |
170 |
GeomCart_old=0.d0 ! Cartesian coordinates |
171 |
|
172 |
Fini=.FALSE. |
173 |
DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) |
174 |
IOpt=IOpt+1 |
175 |
|
176 |
! Calculate the energy and gradient |
177 |
IF (debug) THEN |
178 |
WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt |
179 |
WRITE(*,*) "Energy and Coordinates:" |
180 |
WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12)) |
181 |
WRITE(*,*) "Geom Old:" |
182 |
WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13)) |
183 |
WRITE(*,*) "Grad Old:" |
184 |
WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13)) |
185 |
END IF |
186 |
|
187 |
!Call EGrad(E,Geom,GradTmp,NCoord) |
188 |
IF (COORD.EQ.'BAKER') THEN |
189 |
! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT). |
190 |
! GeomCart has INTENT(OUT) |
191 |
! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable |
192 |
GeomRef=GeomOld |
193 |
Call EGrad_baker(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, & |
194 |
GeomRef,GeomCart_old) |
195 |
GeomCart_old=GeomCart |
196 |
ELSE |
197 |
Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart) |
198 |
END IF |
199 |
|
200 |
IF (debug) THEN |
201 |
WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt |
202 |
WRITE(*,*) "Energy and Coordinates:" |
203 |
WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12)) |
204 |
WRITE(*,*) "Geom Old:" |
205 |
WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13)) |
206 |
END IF |
207 |
|
208 |
IF (IOpt.GE.2) THEN |
209 |
! We have enough geometries and gradient to update the hessian (or its |
210 |
! inverse) |
211 |
GradTmp2=GradTmp-GradOld |
212 |
GeomTmp2=Geom-GeomOld |
213 |
|
214 |
if (debug) Write(*,*) "Call Hupdate_all, Geom" |
215 |
IF (debug) THEN |
216 |
WRITE(*,*) "dx before calling",SHAPE(GeomTmp2) |
217 |
WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord) |
218 |
WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2) |
219 |
WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord) |
220 |
END IF |
221 |
Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local) |
222 |
END IF ! matches IF (IOpt.GE.2) THEN |
223 |
|
224 |
GradOld=GradTmp |
225 |
GeomOld=Geom |
226 |
|
227 |
! GradTmp becomes Step in Step_RFO_all. |
228 |
SELECT CASE (Step_method) |
229 |
CASE ('RFO') |
230 |
Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec) |
231 |
CASE ('GDIIS') |
232 |
Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST) |
233 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1} |
234 |
Geom=0.d0 |
235 |
DO I=1, NCoord |
236 |
! If Hinv=.False., then we need to invert Hess_local |
237 |
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) |
238 |
END DO |
239 |
Geom(:) = NewGeom(:) - Geom(:) |
240 |
! GradTmp now becomes "step" below: |
241 |
GradTmp = Geom - GeomOld |
242 |
CASE ('GDIIST') |
243 |
Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST) |
244 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1} |
245 |
Geom=0.d0 |
246 |
DO I=1, NCoord |
247 |
! If Hinv=.False., then we need to invert Hess_local |
248 |
Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I) |
249 |
END DO |
250 |
Geom(:) = NewGeom(:) - Geom(:) |
251 |
! GradTmp now becomes "step" below: |
252 |
GradTmp = Geom - GeomOld |
253 |
CASE ('GEDIIS') |
254 |
Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST) |
255 |
! GradTmp is actually "step" below: |
256 |
GradTmp = NewGeom - GeomOld |
257 |
!Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec) |
258 |
CASE DEFAULT |
259 |
WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop" |
260 |
STOP |
261 |
END SELECT |
262 |
|
263 |
Fini=.TRUE. |
264 |
! If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0 |
265 |
NormStep=sqrt(DOT_Product(GradTmp,GradTmp)) |
266 |
WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh |
267 |
FactStep=min(1.d0,SMax/NormStep) |
268 |
Fini=(NormStep.LE.SThresh) |
269 |
NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients. |
270 |
WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh |
271 |
Fini=(Fini.AND.(NormStep.LE.GThresh)) |
272 |
if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep |
273 |
GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all. |
274 |
|
275 |
Call Check_step(IGeom,Coord,Nat,NCoord,GeomOld,GradTmp) |
276 |
|
277 |
Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all. |
278 |
|
279 |
IF (debug) THEN |
280 |
WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt |
281 |
SELECT CASE (COORD) |
282 |
CASE ('ZMAT','BAKER') |
283 |
!WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord)) |
284 |
CASE('CART','HYBRID') |
285 |
DO Iat=1,Nat |
286 |
! PFL 30 Apr 2008 -> |
287 |
! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid. |
288 |
! WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |
289 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), & |
290 |
Geom(3*Iat-2:3*Iat) |
291 |
! <- PFL 30 Apr 2008 |
292 |
END DO |
293 |
CASE('MIXED') |
294 |
DO Iat=1,NCart |
295 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |
296 |
Geom(3*Iat-2:3*Iat) |
297 |
END DO |
298 |
|
299 |
SELECT CASE (NCart) |
300 |
CASE(1) |
301 |
if (NAt.GE.2) THEN |
302 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), & |
303 |
IndZmat(2,2),Geom(4) |
304 |
IBEG=5 |
305 |
END IF |
306 |
IF (NAT.GE.3) THEN |
307 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |
308 |
IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6) |
309 |
IBeg=7 |
310 |
END IF |
311 |
CASE(2) |
312 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), & |
313 |
IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8) |
314 |
IBeg=9 |
315 |
CASE DEFAULT |
316 |
IBeg=3*NCart+1 |
317 |
END SELECT |
318 |
|
319 |
DO Iat=max(4,NCart),Nat |
320 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), & |
321 |
IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), & |
322 |
IndZmat(3,3),Geom(IBeg+2)*180./pi |
323 |
IBeg=IBeg+3 |
324 |
END DO |
325 |
|
326 |
CASE DEFAULT |
327 |
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294." |
328 |
STOP |
329 |
END SELECT |
330 |
END IF ! matches IF (debug) THEN |
331 |
|
332 |
END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini)) |
333 |
|
334 |
DEALLOCATE(GradTmp2,GeomTmp2,GradTmp) |
335 |
DEALLOCATE(GradOld,GeomOld) |
336 |
DEALLOCATE(Hess_local) |
337 |
DEALLOCATE(GeomCart_old) |
338 |
DEALLOCATE(NewGeom,NewGradTmp) |
339 |
DEALLOCATE(Hess_local_inv) |
340 |
|
341 |
IF (Fini) THEN |
342 |
WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E |
343 |
ELSE |
344 |
WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom,"*NOT* converged in ",Iopt," cycles, Last Energy = ",E |
345 |
END IF |
346 |
|
347 |
WRITE(*,*) "Initial Geometry:" |
348 |
DO I=1, Nat |
349 |
WRITE(*,'(3(1X,F8.4))') XyzGeomI(IGeom,:,I) |
350 |
END DO |
351 |
WRITE(*,*) "Final Geometry:" |
352 |
DO I=1, Nat |
353 |
WRITE(*,'(3(1X,F8.4))') GeomCart(I,:) |
354 |
!IF (I .GT. 1) Then |
355 |
! WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))& |
356 |
! + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))& |
357 |
! + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3)))) |
358 |
!END IF |
359 |
END DO |
360 |
|
361 |
IF (COORD .EQ. "BAKER") Then |
362 |
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
363 |
ScanCoord=>Coordinate |
364 |
DO WHILE (Associated(ScanCoord%next)) |
365 |
I=I+1 |
366 |
SELECT CASE (ScanCoord%Type) |
367 |
CASE ('BOND') |
368 |
WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I) |
369 |
CASE ('ANGLE') |
370 |
WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, & |
371 |
ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi |
372 |
CASE ('DIHEDRAL') |
373 |
WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,& |
374 |
ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi |
375 |
END SELECT |
376 |
ScanCoord => ScanCoord%next |
377 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
378 |
END IF |
379 |
|
380 |
DEALLOCATE(GeomCart) |
381 |
if (debug) WRITE(*,*) "================== Geom Optimization Over=================" |
382 |
END SUBROUTINE Opt_geom |