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