1      SUBROUTINE DIIS(XP, XPARAM, GP, GRAD, HP, HEAT, HS, NVAR, FRST)
2      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3      INCLUDE 'SIZES'
4      DIMENSION XP(NVAR), XPARAM(NVAR), GP(NVAR),
5     1GRAD(NVAR), HS(NVAR*NVAR)
6      LOGICAL FRST
7************************************************************************
8*                                                                      *
9*     DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE         *
10*                                                                      *
11*     THIS INVOLVES SOLVING FOR C IN XPARAM(NEW) = XPARAM' - HG'       *
12*                                                                      *
13*  WHERE XPARAM' = SUM(C(I)XPARAM(I), THE C COEFFICIENTES COMING FROM  *
14*                                                                      *
15*                   | B   1 | . | C | = | 0 |                          *
16*                   | 1   0 |   |-L |   | 1 |                          *
17*                                                                      *
18* WHERE B(I,J) =GRAD(I)H(T)HGRAD(J)  GRAD(I) = GRADIENT ON CYCLE I     *
19*                              H    = INVERSE HESSIAN                  *
20*                                                                      *
21*                          REFERENCE                                   *
22*                                                                      *
23*  P. CSASZAR, P. PULAY, J. MOL. STRUCT. (THEOCHEM), 114, 31 (1984)    *
24*                                                                      *
25************************************************************************
26************************************************************************
27*                                                                      *
28*     GEOMETRY OPTIMIZATION USING THE METHOD OF DIRECT INVERSION IN    *
29*     THE ITERATIVE SUBSPACE (GDIIS), COMBINED WITH THE BFGS OPTIMIZER *
30*     (A VARIABLE METRIC METHOD)                                       *
31*                                                                      *
32*     WRITTEN BY PETER L. CUMMINS, UNIVERSITY OF SYDNEY, AUSTRALIA     *
33*                                                                      *
34*                              REFERENCE                               *
35*                                                                      *
36*      "COMPUTATIONAL STRATEGIES FOR THE OPTIMIZATION OF EQUILIBRIUM   *
37*     GEOMETRIES AND TRANSITION-STATE STRUCTURES AT THE SEMIEMPIRICAL  *
38*     LEVEL", PETER L. CUMMINS, JILL E. GREADY, J. COMP. CHEM., 10,    *
39*     939-950 (1989).                                                  *
40*                                                                      *
41*     MODIFIED BY JJPS TO CONFORM TO EXISTING MOPAC CONVENTIONS        *
42*                                                                      *
43************************************************************************
44      COMMON /KEYWRD/ KEYWRD
45      PARAMETER (MRESET=15, M2=(MRESET+1)*(MRESET+1))
46      DIMENSION  XSET(MRESET*MAXPAR),GSET(MRESET*MAXPAR), ESET(MRESET)
47      DIMENSION DX(MAXPAR),GSAVE(MAXPAR),
48     1 ERR(MRESET*MAXPAR),B(M2),BS(M2),BST(M2)
49      LOGICAL DEBUG, PRINT
50      CHARACTER*241 KEYWRD
51      DEBUG=.FALSE.
52      PRINT=(INDEX(KEYWRD,' DIIS').NE.0)
53      IF (PRINT) DEBUG=(INDEX(KEYWRD,'DEBUG').NE.0)
54      IF (PRINT)  WRITE(6,'(/,''      ***** BEGIN GDIIS ***** '')')
55C
56C  SPACE SIMPLY LOADS THE CURRENT VALUES OF XPARAM AND GNORM INTO
57C  THE ARRAYS XSET AND GSET
58C
59      CALL SPACE(MRESET,MSET,XPARAM, GRAD, HEAT, NVAR, XSET, GSET, ESET
60     1, FRST)
61C
62C     INITIALIZE SOME VARIABLES AND CONSTANTS
63C
64      NDIIS = MSET
65      MPLUS = MSET + 1
66      MM = MPLUS * MPLUS
67C
68C     COMPUTE THE APPROXIMATE ERROR VECTORS
69C
70      INV=-NVAR
71      DO 30 I=1,MSET
72         INV = INV + NVAR
73         DO 30 J=1,NVAR
74            S = 0.D0
75            KJ=(J*(J-1))/2
76            DO 10 K=1,J
77               KJ = KJ+1
78   10       S = S - HS(KJ) * GSET(INV+K)
79            DO 20 K=J+1,NVAR
80               KJ = (K*(K-1))/2+J
81   20       S = S - HS(KJ) * GSET(INV+K)
82   30 ERR(INV+J) = S
83C
84C     CONSTRUCT THE GDIIS MATRIX
85C
86      DO 40 I=1,MM
87   40 B(I) = 1.D0
88      JJ=0
89      INV=-NVAR
90      DO 50 I=1,MSET
91         INV=INV+NVAR
92         JNV=-NVAR
93         DO 50 J=1,MSET
94            JNV=JNV+NVAR
95            JJ = JJ + 1
96            B(JJ)=0.D0
97            DO 50 K=1,NVAR
98   50 B(JJ) = B(JJ) + ERR(INV+K) * ERR(JNV+K)
99C
100      DO 60 I=MSET-1,1,-1
101         DO 60 J=MSET,1,-1
102   60 B(I*MSET+J+I) = B(I*MSET+J)
103      DO 70 I=1,MPLUS
104         B(MPLUS*I) = 1.D0
105   70 B(MPLUS*MSET+I) = 1.D0
106      B(MM) = 0.D0
107C
108C     ELIMINATE ERROR VECTORS WITH THE LARGEST NORM
109C
110   80 CONTINUE
111      DO 90 I=1,MM
112   90 BS(I) = B(I)
113      IF (NDIIS .EQ. MSET) GO TO 140
114      DO 130 II=1,MSET-NDIIS
115         XMAX = -1.D10
116         ITERA = 0
117         DO 110 I=1,MSET
118            XNORM = 0.D0
119            INV = (I-1) * MPLUS
120            DO 100 J=1,MSET
121  100       XNORM = XNORM + ABS(B(INV + J))
122            IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN
123               XMAX = XNORM
124               ITERA = I
125               IONE = INV + I
126            ENDIF
127  110    CONTINUE
128         DO 120 I=1,MPLUS
129            INV = (I-1) * MPLUS
130            DO 120 J=1,MPLUS
131               JNV = (J-1) * MPLUS
132               IF (J.EQ.ITERA) B(INV + J) = 0.D0
133               B(JNV + I) = B(INV + J)
134  120    CONTINUE
135         B(IONE) = 1.0D0
136  130 CONTINUE
137  140 CONTINUE
138C
139      IF (DEBUG) THEN
140C
141C     OUTPUT THE GDIIS MATRIX
142C
143         WRITE(*,'(/5X,'' GDIIS MATRIX'')')
144         IJ = 0
145         DO 150 I=1,MPLUS
146            DO 150 J=1,I
147               IJ = IJ + 1
148  150    BST(IJ) = B( MPLUS * (J-1) + I)
149         CALL VECPRT(BST,MPLUS)
150      ENDIF
151C
152C     SCALE DIIS MATRIX BEFORE INVERSION
153C
154      DO 160 I=1,MPLUS
155         II = MPLUS * (I-1) + I
156  160 GSAVE(I) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))
157      GSAVE(MPLUS) = 1.D0
158      DO 170 I=1,MPLUS
159         DO 170 J=1,MPLUS
160            IJ = MPLUS * (I-1) + J
161  170 B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J)
162C
163      IF (DEBUG) THEN
164C
165C     OUTPUT SCALED GDIIS MATRIX
166C
167         WRITE(*,'(/5X,'' GDIIS MATRIX (SCALED)'')')
168         IJ = 0
169         DO 180 I=1,MPLUS
170            DO 180 J=1,I
171               IJ = IJ + 1
172  180    BST(IJ) = B( MPLUS * (J-1) + I)
173         CALL VECPRT(BST,MPLUS)
174      ENDIF
175C
176C     INVERT THE GDIIS MATRIX
177C
178      CALL MINV(B,MPLUS,DET)
179C
180      DO 190 I=1,MPLUS
181         DO 190 J=1,MPLUS
182            IJ = MPLUS * (I-1) + J
183  190 B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J)
184C
185C     COMPUTE THE INTERMEDIATE INTERPOLATED PARAMETER AND GRADIENT
186C     VECTORS
187C
188      DO 200 K=1,NVAR
189         XP(K) = 0.D0
190         GP(K) = 0.D0
191         DO 200 I=1,MSET
192            INK = (I-1) * NVAR + K
193            XP(K) = XP(K) + B(MPLUS*MSET+I) * XSET(INK)
194  200 GP(K) = GP(K) + B(MPLUS*MSET+I) * GSET(INK)
195      HP=0.D0
196      DO 210 I=1,MSET
197  210 HP=HP+B(MPLUS*MSET+I)*ESET(I)
198C
199      DO 220 K=1,NVAR
200  220 DX(K) = XPARAM(K) - XP(K)
201      XNORM = SQRT(DOT(DX,DX,NVAR))
202      IF (PRINT) THEN
203         WRITE (6,'(/10X,''DEVIATION IN X '',F7.4,8X,''DETERMINANT '',
204     1 G9.3)') XNORM,DET
205         WRITE(6,'(10X,''GDIIS COEFFICIENTS'')')
206         WRITE(6,'(10X,5F12.5)') (B(MPLUS*MSET+I),I=1,MSET)
207      ENDIF
208C
209C     THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY!
210C
211      THRES = MAX(10.D0**(-NVAR), 1.D-25)
212      IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN
213         IF (PRINT) WRITE(6,'(10X,''THE DIIS MATRIX IS ILL CONDITIONED''
214     1, /10X,'' - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - '',
215     2 /10X,''THE DIIS STEP WILL BE REPEATED WITH A SMALLER SPACE'')')
216         DO 230 K=1,MM
217  230    B(K) = BS(K)
218         NDIIS = NDIIS - 1
219         IF (NDIIS .GT. 0) GO TO 80
220         IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')')
221         DO 240 K=1,NVAR
222            XP(K) = XPARAM(K)
223  240    GP(K) = GRAD(K)
224C
225      ENDIF
226      IF (PRINT)  WRITE(6,'(/,''      ***** END GDIIS ***** '',/)')
227C
228      RETURN
229      END
230      SUBROUTINE SPACE(MRESET, MSET, XPARAM, GRAD, HEAT, NVAR,
231     1XSET, GSET, ESET, FRST)
232      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
233      INCLUDE 'SIZES'
234      DIMENSION XPARAM(NVAR), GRAD(NVAR)
235      DIMENSION XSET(MRESET*NVAR),GSET(MRESET*NVAR), ESET(MRESET)
236      LOGICAL FRST
237C
238C     UPDATE PARAMETER AND GRADIENT SUBSPACE
239C
240      IF(FRST)THEN
241         NRESET=MIN(NVAR/2,MRESET)
242         FRST=.FALSE.
243         MSET=0
244      ENDIF
245C
246      IF (MSET .EQ. NRESET) THEN
247         DO 10 I=1,MSET-1
248            MI = NVAR*(I-1)
249            NI = NVAR*I
250            ESET(I)=ESET(I+1)
251            DO 10 K=1,NVAR
252               XSET(MI+K) = XSET(NI+K)
253   10    GSET(MI+K) = GSET(NI+K)
254         MSET=NRESET-1
255      ENDIF
256C
257C     STORE THE CURRENT POINT
258C
259      DO 20 K=1,NVAR
260         NMK = NVAR*MSET+K
261         XSET(NMK) = XPARAM(K)
262   20 GSET(NMK) = GRAD(K)
263      MSET=MSET+1
264      ESET(MSET)=HEAT
265C
266      RETURN
267      END
268      SUBROUTINE MINV(A,N,D)
269      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
270      INCLUDE 'SIZES'
271      DIMENSION A(*)
272**********************************************************************
273*
274*     INVERT A MATRIX USING GAUSS-JORDAN METHOD.  PART OF DIIS
275*     A - INPUT MATRIX (MUST BE A GENERAL MATRIX), DESTROYED IN
276*        COMPUTATION AND REPLACED BY RESULTANT INVERSE.
277*     N - ORDER OF MATRIX A
278*     D - RESULTANT DETERMINANT
279*
280**********************************************************************
281      DIMENSION M(MAXPAR), L(MAXPAR)
282C
283C     SEARCH FOR LARGEST ELEMENT
284C
285      D=1.0D0
286      NK=-N
287      DO 180 K=1,N
288         NK=NK+N
289         L(K)=K
290         M(K)=K
291         KK=NK+K
292         BIGA=A(KK)
293         DO 20 J=K,N
294            IZ=N*(J-1)
295            DO 20 I=K,N
296               IJ=IZ+I
297   10          IF (ABS(BIGA).LT.ABS(A(IJ)))THEN
298                  BIGA=A(IJ)
299                  L(K)=I
300                  M(K)=J
301               ENDIF
302   20    CONTINUE
303C
304C     INTERCHANGE ROWS
305C
306         J=L(K)
307         IF (J-K) 50,50,30
308   30    KI=K-N
309         DO 40 I=1,N
310            KI=KI+N
311            HOLD=-A(KI)
312            JI=KI-K+J
313            A(KI)=A(JI)
314   40    A(JI)=HOLD
315C
316C     INTERCHANGE COLUMNS
317C
318   50    I=M(K)
319         IF (I-K) 80,80,60
320   60    JP=N*(I-1)
321         DO 70 J=1,N
322            JK=NK+J
323            JI=JP+J
324            HOLD=-A(JK)
325            A(JK)=A(JI)
326   70    A(JI)=HOLD
327C
328C     DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS
329C     CONTAINED IN BIGA)
330C
331   80    IF (BIGA) 100,90,100
332   90    D=0.0D0
333         RETURN
334  100    DO 120 I=1,N
335            IF (I-K) 110,120,110
336  110       IK=NK+I
337            A(IK)=A(IK)/(-BIGA)
338  120    CONTINUE
339C  REDUCE MATRIX
340         DO 150 I=1,N
341            IK=NK+I
342            HOLD=A(IK)
343            IJ=I-N
344            DO 150 J=1,N
345               IJ=IJ+N
346               IF (I-K) 130,150,130
347  130          IF (J-K) 140,150,140
348  140          KJ=IJ-I+K
349               A(IJ)=HOLD*A(KJ)+A(IJ)
350  150    CONTINUE
351C
352C     DIVIDE ROW BY PIVOT
353C
354         KJ=K-N
355         DO 170 J=1,N
356            KJ=KJ+N
357            IF (J-K) 160,170,160
358  160       A(KJ)=A(KJ)/BIGA
359  170    CONTINUE
360C
361C     PRODUCT OF PIVOTS
362C
363         D=MAX(-1.D25,MIN(1.D25,D))
364         D=D*BIGA
365C
366C     REPLACE PIVOT BY RECIPROCAL
367C
368         A(KK)=1.0D0/BIGA
369  180 CONTINUE
370C
371C     FINAL ROW AND COLUMN INTERCHANGE
372C
373      K=N
374  190 K=(K-1)
375      IF (K) 260,260,200
376  200 I=L(K)
377      IF (I-K) 230,230,210
378  210 JQ=N*(K-1)
379      JR=N*(I-1)
380      DO 220 J=1,N
381         JK=JQ+J
382         HOLD=A(JK)
383         JI=JR+J
384         A(JK)=-A(JI)
385  220 A(JI)=HOLD
386  230 J=M(K)
387      IF (J-K) 190,190,240
388  240 KI=K-N
389      DO 250 I=1,N
390         KI=KI+N
391         HOLD=A(KI)
392         JI=KI-K+J
393         A(KI)=-A(JI)
394  250 A(JI) =HOLD
395      GO TO 190
396  260 RETURN
397      END
398