root / src / Step_DIIS_all.f90 @ 8
History  View  Annotate  Download (16.3 kB)
1 
!C HEAT is never used, not even in call of Space(...) 

2 
!C Geom = input parameter vector (Geometry). 
3 
!C Grad = input gradient vector. 
4 
!C Geom_new = New Geometry. 
5  
6 
SUBROUTINE Step_diis_all(NGeomF,IGeom,Step,Geom,Grad,HP,HEAT,Hess,NCoord,allocation_flag,Tangent) 
7 
! IMPLICIT DOUBLE PRECISION (AH,OZ) 
8  
9 
USE Io_module, only : IoOut 
10 
USE Path_module, only : Vfree 
11 

12 
IMPLICIT NONE 
13 
INTEGER, parameter :: KINT = kind(1) 
14 
INTEGER, parameter :: KREAL = kind(1.0d0) 
15 

16 
! INCLUDE 'SIZES' 
17 

18 
INTEGER(KINT) :: NGeomF,IGeom 
19 
INTEGER(KINT), INTENT(IN) :: NCoord 
20 

21 
REAL(KREAL) :: Geom(NCoord), Grad(NCoord) 
22 
REAL(KREAL) :: Hess(NCoord*NCoord),Step(NCoord) 
23 
REAL(KREAL) :: HEAT,HP 
24 
LOGICAL :: allocation_flag 
25 
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) 
26  
27 
!************************************************************************ 
28 
!* * 
29 
!* DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE * 
30 
!* * 
31 
!* THIS INVOLVES SOLVING FOR C IN Geom(NEW) = Geom'  HG' * 
32 
!* * 
33 
!* WHERE Geom' = SUM(C(I)Geom(I), THE C COEFFICIENTES COMING FROM * 
34 
!* * 
35 
!*  B 1  .  C  =  0  * 
36 
!*  1 0  L   1  * 
37 
!* * 
38 
!* WHERE B(I,J) =GRAD(I)H(T)HGRAD(J) GRAD(I) = GRADIENT ON CYCLE I * 
39 
!* Hess = INVERSE HESSIAN * 
40 
!* * 
41 
!* REFERENCE * 
42 
!* * 
43 
!* P. CSASZAR, P. PULAY, J. MOL. STRUCT. (THEOCHEM), 114, 31 (1984) * 
44 
!* * 
45 
!************************************************************************ 
46 
!************************************************************************ 
47 
!* * 
48 
!* GEOMETRY OPTIMIZATION USING THE METHOD OF DIRECT INVERSION IN * 
49 
!* THE ITERATIVE SUBSPACE (GDIIS), COMBINED WITH THE BFGS OPTIMIZER * 
50 
!* (A VARIABLE METRIC METHOD) * 
51 
!* * 
52 
!* WRITTEN BY PETER L. CUMMINS, UNIVERSITY OF SYDNEY, AUSTRALIA * 
53 
!* * 
54 
!* REFERENCE * 
55 
!* * 
56 
!* "COMPUTATIONAL STRATEGIES FOR THE OPTIMIZATION OF EQUILIBRIUM * 
57 
!* GEOMETRIES AND TRANSITIONSTATE STRUCTURES AT THE SEMIEMPIRICAL * 
58 
!* LEVEL", PETER L. CUMMINS, JILL E. GREADY, J. COMP. CHEM., 10, * 
59 
!* 939950 (1989). * 
60 
!* * 
61 
!* MODIFIED BY JJPS TO CONFORM TO EXISTING MOPAC CONVENTIONS * 
62 
!* * 
63 
!************************************************************************ 
64  
65 
! MRESET = maximum number of iterations. 
66 
INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 
67 
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:),GradSet(:,:),ERR(:,:) ! MRESET*NCoord 
68 
REAL(KREAL), SAVE :: ESET(MRESET) 
69 
REAL(KREAL), ALLOCATABLE, SAVE :: DXTMP(:,:),GSAVE(:,:) !NCoord, why DXTMP has SAVE attribute?? 
70 
REAL(KREAL) :: B(M2), BS(M2) 
71 
LOGICAL DEBUG, PRINT 
72 
INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) 
73 
LOGICAL, ALLOCATABLE, SAVE :: FRST(:) 
74 
INTEGER(KINT) :: NDIIS, MPLUS, INV, ITERA, MM, NFree, I, J, K 
75 
INTEGER(KINT) :: JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, Isch, Idx 
76 
REAL(KREAL) :: XMax, XNorm, S, DET, THRES, Norm 
77 
REAL(KREAL), PARAMETER :: eps=1e12 
78 
REAL(KREAL), PARAMETER :: crit=1e8 
79 
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord 
80 
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree 
81 
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree 
82 
REAL(KREAL), ALLOCATABLE :: Grad_free(:),Grad_new_free_inter(:),Step_free(:) ! NFree 
83 
REAL(KREAL), ALLOCATABLE :: Geom_free(:),Geom_new_free_inter(:),Geom_new_free(:) ! NFree 
84 
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:),GradSet_free(:,:) 
85 

86 
DEBUG=.TRUE. 
87 
PRINT=.TRUE. 
88  
89 
IF (PRINT) WRITE(*,'(/,'' BEGIN GDIIS '')') 
90 

91 
! Initialization 
92 
IF (allocation_flag) THEN ! allocation_flag = .TRUE. at the begining and effective for all geometries in path. 
93 
! FRST(IGeom) will be set to False in Space, so no need to modify it here 
94 
IF (ALLOCATED(GeomSet)) THEN 
95 
IF (PRINT) WRITE(*,'(/,'' In FRST, GDIIS Dealloc '')') 
96 
DEALLOCATE(GeomSet,GradSet,ERR,DXTMP,GSave,GeomSet_free,GradSet_free) 
97 
RETURN 
98 
ELSE 
99 
! these allocated arrays need to be properly deallocated. 
100 
IF (PRINT) WRITE(*,'(/,'' In FRST, GDIIS Alloc '')') 
101 
ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),ERR(NGeomF,MRESET*NCoord)) 
102 
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord)) 
103 
ALLOCATE(DXTMP(NGeomF,NCoord),GSAVE(NGeomF,NCoord),MSET(NGeomF),FRST(NGeomF)) 
104 
DO I=1,NGeomF 
105 
FRST(I) = .TRUE. 
106 
GeomSet(I,:) = 0.d0 
107 
GradSet(I,:) = 0.d0 
108 
ERR(I,:) = 0.d0 
109 
GeomSet_free(I,:) = 0.d0 
110 
GradSet_free(I,:) = 0.d0 
111 
DXTMP(I,:)=0.d0 
112 
GSAVE(I,:)=0.d0 
113 
END DO 
114 
MSET(:)=0 
115 
END IF 
116 
allocation_flag = .FALSE. 
117 
END IF 
118 

119 
! Addded from here: 
120 
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord) 
121 
! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0 
122 
Norm=sqrt(dot_product(Tangent,Tangent)) 
123 
IF (Norm.GT.eps) THEN 
124 
ALLOCATE(Tanf(NCoord)) 
125  
126 
! We normalize Tangent 
127 
Tangent=Tangent/Norm 
128  
129 
! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23) 
130 
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints. 
131 
DO I=1,NCoord 
132 
Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent) 
133 
END DO 
134 
Tangent=0.d0 
135 
DO I=1,NCoord 
136 
Tangent=Tangent+Tanf(I)*Vfree(:,I) 
137 
END DO 
138 
! first we subtract Tangent from vfree 
139 
DO I=1,NCoord 
140 
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) 
141 
Vfree(:,I)=Vfree(:,I)Norm*Tangent 
142 
END DO 
143  
144 
Idx=0 
145 
! Schmidt orthogonalization of the Vfree vectors 
146 
DO I=1,NCoord 
147 
! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable. 
148 
DO Isch=1,2 
149 
DO J=1,Idx 
150 
Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/))) 
151 
Vfree(:,I)=Vfree(:,I)Norm*Vfree(:,J) 
152 
END DO 
153 
END DO 
154 
Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/))) 
155 
IF (Norm.GE.crit) THEN 
156 
Idx=Idx+1 
157 
Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm) 
158 
END IF 
159 
END DO 
160 

161 
Print *, 'Idx=', Idx 
162 
IF (Idx/= NCoord1) THEN 
163 
WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom 
164 
WRITE(IoOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom 
165 
STOP 
166 
END IF 
167 

168 
DEALLOCATE(Tanf) 
169 
NFree=Idx 
170 
ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN 
171 
if (debug) WRITE(*,*) "Tangent=0, using full displacement" 
172 
NFree=NCoord 
173 
END IF !IF (Norm.GT.eps) THEN 
174 

175 
if (debug) WRITE(*,*) 'DBG Step_DIIS_All, IGeom, NFree=', IGeom, NFree 
176  
177 
! We now calculate the new step 
178 
! we project the hessian onto the free vectors 
179 
ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree),Grad_new_free_inter(NFree)) 
180 
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free_inter(NFree),Geom_new_free(NFree)) 
181 
DO J=1,NFree 
182 
DO I=1,NCoord 
183 
Htmp(I,J)=0.d0 
184 
DO K=1,NCoord 
185 
Htmp(I,J)=Htmp(I,J)+Hess(((I1)*NCoord)+K)*Vfree(K,J) 
186 
END DO 
187 
END DO 
188 
END DO 
189 
DO J=1,NFree 
190 
DO I=1,NFree 
191 
HFree(I+((J1)*NFree))=0.d0 
192 
DO K=1,NCoord 
193 
HFree(I+((J1)*NFree))=HFree(I+((J1)*NFree))+Vfree(K,I)*Htmp(K,J) 
194 
END DO 
195 
END DO 
196 
END DO 
197  
198 
DO I=1,NFree 
199 
Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad) 
200 
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom) 
201 
END DO 
202 
!Added Ends here.*********************************************** 
203 

204 
!C SPACE SIMPLY LOADS THE CURRENT VALUES OF Geom AND GRAD INTO 
205 
!C THE ARRAYS GeomSet AND GradSet 
206 
!C HEAT is never used, not even in Space_all(...) 
207  
208 
CALL Space_all(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) 
209 

210 
IF (PRINT) WRITE(*,'(/,'' GDIIS after Space '')') 
211 

212 
DO J=1,MSet(IGeom) 
213 
DO K=1,NFree 
214 
GradSet_free(IGeom,((J1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& 
215 
GradSet(IGeom,((J1)*NCoord)+1:((J1)*NCoord)+NCoord)) 
216 
GeomSet_free(IGeom,((J1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& 
217 
GeomSet(IGeom,((J1)*NCoord)+1:((J1)*NCoord)+NCoord)) 
218 
END DO 
219 
END DO 
220 
!C 
221 
!C INITIALIZE SOME VARIABLES AND CONSTANTS 
222 
!C 
223 
NDIIS = MSET(IGeom) 
224 
MPLUS = MSET(IGeom) + 1 
225 
MM = MPLUS * MPLUS 
226 
!C 
227 
!C COMPUTE THE APPROXIMATE ERROR VECTORS 
228 
!C 
229 
INV=NFree 
230 
DO 30 I=1,MSET(IGeom) 
231 
INV = INV + NFree 
232 
DO 30 J=1,NFree 
233 
S = 0.D0 
234 
KJ=(J*(J1))/2 
235 
DO 10 K=1,J 
236 
KJ = KJ+1 
237 
10 S = S  HFree(KJ) * GradSet_free(IGeom,INV+K) 
238 
DO 20 K=J+1,NFree 
239 
KJ = (K*(K1))/2+J 
240 
20 S = S  HFree(KJ) * GradSet_free(IGeom,INV+K) 
241 
30 ERR(IGeom,INV+J) = S 
242 

243 
!C 
244 
!C CONSTRUCT THE GDIIS MATRIX 
245 
!C 
246 
DO 40 I=1,MM 
247 
40 B(I) = 1.D0 
248 

249 
JJ=0 
250 
INV=NFree 
251 
DO 50 I=1,MSET(IGeom) 
252 
INV=INV+NFree 
253 
JNV=NFree 
254 
DO 50 J=1,MSET(IGeom) 
255 
JNV=JNV+NFree 
256 
JJ = JJ + 1 
257 
B(JJ)=0.D0 
258 
DO 50 K=1,NFree 
259 
!Print *, 'B(',JJ,')=', B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K) 
260 
50 B(JJ) = B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K) 
261  
262 
! The following shifting is required to correct indices of B_ij elements in the GDIIS matrix. 
263 
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero. 
264 
DO 60 I=MSET(IGeom)1,1,1 
265 
DO 60 J=MSET(IGeom),1,1 
266 
60 B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) 
267 

268 
! For the last row and last column of GEDIIS matrix: 
269 
DO 70 I=1,MPLUS 
270 
B(MPLUS*I) = 1.D0 
271 
70 B(MPLUS*MSET(IGeom)+I) = 1.D0 
272 
B(MM) = 0.D0 
273 
!C 
274 
!C ELIMINATE ERROR VECTORS WITH THE LARGEST NORM 
275 
!C 
276 
80 CONTINUE 
277 
DO 90 I=1,MM 
278 
90 BS(I) = B(I) 
279 

280 
IF (NDIIS .EQ. MSET(IGeom)) GO TO 140 
281 
DO 130 II=1,MSET(IGeom)NDIIS 
282 
XMAX = 1.D10 
283 
ITERA = 0 
284 
DO 110 I=1,MSET(IGeom) 
285 
XNORM = 0.D0 
286 
INV = (I1) * MPLUS 
287 
DO 100 J=1,MSET(IGeom) 
288 
100 XNORM = XNORM + ABS(B(INV + J)) 
289 
IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN 
290 
XMAX = XNORM 
291 
ITERA = I 
292 
IONE = INV + I 
293 
ENDIF 
294 
110 CONTINUE 
295 

296 
DO 120 I=1,MPLUS 
297 
INV = (I1) * MPLUS 
298 
DO 120 J=1,MPLUS 
299 
JNV = (J1) * MPLUS 
300 
IF (J.EQ.ITERA) B(INV + J) = 0.D0 
301 
B(JNV + I) = B(INV + J) 
302 
!Print *,'B(JNV + I)=',B(JNV + I) 
303 
120 CONTINUE 
304 
B(IONE) = 1.0D0 
305 
130 CONTINUE 
306 
140 CONTINUE 
307 
!C 
308 
!C OUTPUT THE GDIIS MATRIX 
309 
!C 
310 
IF (DEBUG) THEN 
311 
WRITE(*,'(/5X,'' GDIIS MATRIX'')') 
312 
ITmp=min(12,MPLUS) 
313 
DO IJ=1,MPLUS 
314 
WRITE(*,'(12(F12.4,1X))') B((IJ1)*MPLUS+1:(IJ1)*MPLUS+ITmp) 
315 
END DO 
316 
ENDIF 
317 
!C 
318 
!C SCALE DIIS MATRIX BEFORE INVERSION 
319 
!C 
320 
DO 160 I=1,MPLUS 
321 
II = MPLUS * (I1) + I 
322 
!Print *, 'B(',II,')=', B(II) 
323 
!Print *, 'GSave(',IGeom,',',I,')=', 1.D0 / DSQRT(1.D20+DABS(B(II))) 
324 
160 GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D20+DABS(B(II))) 
325 

326 
GSAVE(IGeom,MPLUS) = 1.D0 
327 
!Print *, 'GSave(',IGeom,',',MPlus,')=1.D0' 
328 

329 
DO 170 I=1,MPLUS 
330 
DO 170 J=1,MPLUS 
331 
IJ = MPLUS * (I1) + J 
332 
170 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) 
333 
!C 
334 
!C OUTPUT SCALED GDIIS MATRIX 
335 
!C 
336 
IF (DEBUG) THEN 
337 
WRITE(*,'(/5X,'' GDIIS MATRIX (SCALED)'')') 
338 
ITmp=min(12,MPLUS) 
339 
DO IJ=1,MPLUS 
340 
WRITE(*,'(12(F12.4,1X))') B((IJ1)*MPLUS+1:(IJ1)*MPLUS+ITmp) 
341 
END DO 
342 
ENDIF 
343 
!C 
344 
!C INVERT THE GDIIS MATRIX B 
345 
!C 
346 
CALL MINV(B,MPLUS,DET) ! matrix inversion. 
347  
348 
DO 190 I=1,MPLUS 
349 
DO 190 J=1,MPLUS 
350 
IJ = MPLUS * (I1) + J 
351 
!Print *, 'B(',IJ,')=', B(IJ) 
352 
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I) 
353 
!Print *, 'GSAVE(',IGeom,',',J,')=', GSAVE(IGeom,J) 
354 
!Print *, 'B(',IJ,')=', B(IJ) * GSAVE(I) * GSAVE(J) 
355 
190 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) 
356 
!C 
357 
!C COMPUTE THE INTERMEDIATE INTERPOLATED PARAMETER AND GRADIENT VECTORS 
358 
!C 
359 
!Print *, 'MSET(',IGeom,')=', MSET(IGeom), ' MPLUS=', MPLUS 
360 
DO 200 K=1,NFree 
361 
Geom_new_free_inter(K) = 0.D0 
362 
Grad_new_free_inter(K) = 0.D0 
363 
DO 200 I=1,MSET(IGeom) 
364 
INK = (I1) * NFree + K 
365 
Geom_new_free_inter(K) = Geom_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GeomSet_free(IGeom,INK) 
366 
!Print *, 'Geom_new_free_inter(',K,')=', Geom_new_free_inter(K) 
367 
!Print *, 'B(MPLUS*MSET(',IGeom,')+',I,')=', B(MPLUS*MSET(IGeom)+I) 
368 
!Print *, 'GeomSet_free(',IGeom,',',INK,')=', GeomSet_free(IGeom,INK) 
369 
200 Grad_new_free_inter(K) = Grad_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GradSet_free(IGeom,INK) 
370 
HP=0.D0 
371 
DO 210 I=1,MSET(IGeom) 
372 
210 HP=HP+B(MPLUS*MSET(IGeom)+I)*ESET(I) 
373 
DO 220 K=1,NFree 
374 
220 DXTMP(IGeom,K) = Geom_free(K)  Geom_new_free_inter(K) 
375 
XNORM = SQRT(DOT_PRODUCT(DXTMP(IGeom,1:NFree),DXTMP(IGeom,1:NFree))) 
376 
IF (PRINT) THEN 
377 
WRITE (6,'(/10X,''DEVIATION IN X '',F10.6, 8X,''DETERMINANT '',G9.3)') XNORM,DET 
378 
WRITE(*,'(10X,''GDIIS COEFFICIENTS'')') 
379 
WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom)) 
380 
ENDIF 
381  
382 
!C THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY! 
383 
THRES = MAX(10.D0**(NFree), 1.D25) 
384 
IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN 
385 
IF (PRINT)THEN 
386 
WRITE(*,*) "THE DIIS MATRIX IS ILL CONDITIONED" 
387 
WRITE(*,*) "  PROBABLY, VECTORS ARE LINEARLY DEPENDENT  " 
388 
WRITE(*,*) "THE DIIS STEP WILL BE REPEATED WITH A SMALLER SPACE" 
389 
END IF 
390 
DO 230 K=1,MM 
391 
230 B(K) = BS(K) 
392 
NDIIS = NDIIS  1 
393 
IF (NDIIS .GT. 0) GO TO 80 
394 
IF (PRINT) WRITE(*,'(10X,''NEWTONRAPHSON STEP TAKEN'')') 
395 
DO 240 K=1,NFree 
396 
Geom_new_free_inter(K) = Geom_free(K) 
397 
240 Grad_new_free_inter(K) = Grad_free(K) 
398 
ENDIF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN, L378 
399 

400 
! q_{m+1} = q'_{m+1}  H^{1}g'_{m+1} 
401 
Geom_new_free=0.d0 
402 
DO I = 1, NFree 
403 
DO J = 1, NFree 
404 
! If Hinv=.False., then we need to invert Hess 
405 
!Geom_new_free(:) = Geom_new_free(:) + HFree(:,I)*Grad_new_free_inter(I) 
406 
Geom_new_free(J) = Geom_new_free(J) + HFree(I+((J1)*NFree))*Grad_new_free_inter(I) 
407 
END DO 
408 
END DO 
409 
Geom_new_free(:) = Geom_new_free_inter(:)  Geom_new_free(:) 
410 

411 
Step_free = Geom_new_free  Geom_free 
412 

413 
Step = 0.d0 
414 
DO I=1,NFree 
415 
Step = Step + Step_free(I)*Vfree(:,I) 
416 
END DO 
417  
418 
DEALLOCATE(Hfree,Htmp,Grad_free,Grad_new_free_inter,Step_free,Geom_free) 
419 
DEALLOCATE(Geom_new_free_inter,Geom_new_free) 
420 

421 
IF (PRINT) WRITE(*,'(/,'' END GDIIS '',/)') 
422 

423 
END SUBROUTINE Step_diis_all 