1 SUBROUTINE ITER (H, W, WJ, WK, EE, FULSCF,RAND) 2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 3 INCLUDE 'SIZES' 4 DOUBLE PRECISION MECI 5 DIMENSION H(*), W(*), WJ(*), WK(*) 6 COMMON /FOKMAT/ F(MPACK), FB(MPACK) 7 COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK) 8 COMMON /VECTOR/ C(MORB2),EIGS(MAXORB),CBETA(MORB2),EIGB(MAXORB) 9 COMMON /GRADNT/ DUMY(MAXPAR),GNORM 10 COMMON /LAST / LAST 11 COMMON /MESAGE/ IFLEPO,IITER 12 COMMON /ATHEAT/ ATHEAT 13 COMMON /ENUCLR/ ENUCLR 14 COMMON /CITERM/ XI,XJ,XK 15 COMMON /PATH / LATOM,LPARAM,REACT(200) 16 COMMON /NUMCAL/ NUMCAL 17 COMMON /SCFTYP/ EMIN, LIMSCF 18C ***** Modified by Jiro Toyoda at 1994-05-25 ***** 19C COMMON /TIME / TIME0 20 COMMON /TIMEC / TIME0 21C ***************************** at 1994-05-25 ***** 22 LOGICAL FULSCF, RAND, LIMSCF 23 DOUBLE PRECISION WJ, WK 24C*********************************************************************** 25C 26C ITER GENERATES A SCF FIELD AND RETURNS THE ENERGY IN "ENERGY" 27C 28C THE MAIN ARRAYS USED IN ITER ARE: 29C P ONLY EVER CONTAINS THE TOTAL DENSITY MATRIX 30C PA ONLY EVER CONTAINS THE ALPHA DENSITY MATRIX 31C PB ONLY EVER CONTAINS THE BETA DENSITY MATRIX 32C C ONLY EVER CONTAINS THE EIGENVECTORS 33C H ONLY EVER CONTAINS THE ONE-ELECTRON MATRIX 34C F STARTS OFF CONTAINING THE ONE-ELECTRON MATRIX, 35C AND IS USED TO HOLD THE FOCK MATRIX 36C W ONLY EVER CONTAINS THE TWO-ELECTRON MATRIX 37C 38C THE MAIN INTEGERS CONSTANTS IN ITER ARE: 39C 40C LINEAR SIZE OF PACKED TRIANGLE = NORBS*(NORBS+1)/2 41C 42C THE MAIN INTEGER VARIABLES ARE 43C NITER NUMBER OF ITERATIONS EXECUTED 44C 45C PRINCIPAL REFERENCES: 46C 47C ON MNDO: "GROUND STATES OF MOLECULES. 38. THE MNDO METHOD. 48C APPROXIMATIONS AND PARAMETERS." 49C DEWAR, M.J.S., THIEL,W., J. AM. CHEM. SOC.,99,4899,(1977). 50C ON SHIFT: "THE DYNAMIC 'LEVEL SHIFT' METHOD FOR IMPROVING THE 51C CONVERGENCE OF THE SCF PROCEDURE", A. V. MITIN, J. COMP. 52C CHEM. 9, 107-110 (1988) 53C ON HALF-ELECTRON: "MINDO/3 COMPARISON OF THE GENERALIZED S.C.F. 54C COUPLING OPERATOR AND "HALF-ELECTRON" METHODS FOR 55C CALCULATING THE ENERGIES AND GEOMETRIES OF OPEN SHELL 56C SYSTEMS" 57C DEWAR, M.J.S., OLIVELLA, S., J. CHEM. SOC. FARA. II, 58C 75,829,(1979). 59C ON PULAY'S CONVERGER: "CONVERGANCE ACCELERATION OF ITERATIVE 60C SEQUENCES. THE CASE OF SCF ITERATION", PULAY, P., 61C CHEM. PHYS. LETT, 73, 393, (1980). 62C ON CNVG: IT ENCORPORATES THE IMPROVED ITERATION SCHEME (IIS) BY 63C PIOTR BADZIAG & FRITZ SOLMS. ACCEPTED FOR PUBLISHING 64C IN COMPUTERS & CHEMISTRY 65C ON PSEUDODIAGONALISATION: "FAST SEMIEMPIRICAL CALCULATIONS", 66C STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM., 67C 3, 227, (1982) 68C 69C*********************************************************************** 70 DIMENSION POLD(MPACK), POLD2(MPACK), POLD3(MAXORB+400) 71 DIMENSION PBOLD(MPACK), PBOLD2(MPACK), PBOLD3(MAXORB+400) 72************************************************************************ 73* * 74* PACK ALL THE ARRAYS USED BY PULAY INTO A COMMON BLOCK SO THAT THEY * 75* CAN BE USED BY THE C.I. DERIVATIVE, IF NEEDED * 76* * 77************************************************************************ 78 COMMON /WORK3/POLD,POLD2,PBOLD,PBOLD2 79 COMMON /WORK1/ AR1,AR2,AR3,AR4,BR1,BR2,BR3,BR4 80 DIMENSION AR1(2*NPULAY), AR2(2*NPULAY), AR3(2*NPULAY), 81 1 AR4(2*NPULAY) 82 DIMENSION BR1(2*NPULAY), BR2(2*NPULAY), BR3(2*NPULAY), 83 1 BR4(7*NPULAY) 84 DIMENSION ESCF0(10) 85 COMMON /PRECI / SELCON 86 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 87 1 NLAST(NUMATM), NORBS, NELECS, 88 2 NALPHA, NBETA, NCLOSE, NOPEN, NDUMY, FRACT 89 3 /MOLORB/ DUMMY(MAXORB),PDIAG(MAXORB) 90 4 /KEYWRD/ KEYWRD 91 5 /NUMSCF/ NSCF 92 SAVE ICALCN, DEBUG, PRTFOK, PRTEIG, PRTDEN, PRT1EL, ABPRT 93 SAVE LINEAR, MINPRT, NEWDG, SCFCRT, PRTPL, PRTVEC, PL 94 SAVE BSHIFT, PLTEST, ITRMAX, NA2EL, NA1EL, NB2EL,NB1EL 95 SAVE IFILL, CAMKIN, CI, OKPULY, UHF, SCF1, OKNEWD, TIMES 96 SAVE FORCE, ALLCON, TRANS, HALFE, W1, W2, RANDOM 97 CHARACTER KEYWRD*241, ABPRT(3)*5, GETNAM*80 98 LOGICAL PRTFOK,PRTEIG,PRTDEN, DEBUG, TIMES, CI 99 1,UHF, NEWDG, SCF1, HALFE, FORCE, PRT1EL,PRTPL, OKNEWD 100 2,MINPRT, FRST, BFRST, OKPULY, READY, PRTVEC, 101 3CAMKIN, ALLCON, MAKEA, MAKEB, INCITR, CAPPS, TIMITR 102 DATA ICALCN/0/, DEBUG/.FALSE./, PRTFOK/.FALSE./ 103 DATA PRTEIG/.FALSE./,PRTDEN/.FALSE./ 104 DATA PRT1EL/.FALSE./ 105 DATA ABPRT/' ','ALPHA',' BETA'/ 106C 107C INITIALIZE 108C 109 IFILL=0 110 IHOMO=MAX(1,NCLOSE+NALPHA) 111 IHOMOB=MAX(1,NCLOSE+NBETA) 112 EOLD=1.D2 113 READY=.FALSE. 114 IF (ICALCN.NE.NUMCAL) THEN 115 CALL EPSETA(EPS,ETA) 116C 117C ULTIMATE SCF CRITERION: HEAT OF FORMATION CONVERGED WITHIN A FACTOR 118C OF 10 OF THE LIMITING PRECISION OF THE COMPUTER 119C 120 EPS=23.061D0*EPS*10.D0 121 IRRR=5 122 SHIFT=0.D0 123 ICALCN=NUMCAL 124 SHFMAX=20.D0 125 LINEAR=(NORBS*(NORBS+1))/2 126C 127C DEBUG KEY-WORDS WORKED OUT 128C 129 DEBUG=( INDEX(KEYWRD,'DEBUG') .NE. 0 ) 130 MINPRT=(INDEX(KEYWRD,'SADDLE')+ 131 1 LATOM .EQ.0 .OR. DEBUG) 132 PRTEIG=( INDEX(KEYWRD,'EIGS') .NE. 0 ) 133 PRTPL =( INDEX(KEYWRD,' PL ') .NE.0 ) 134 PRT1EL=( INDEX(KEYWRD,'1ELE') .NE.0 .AND. DEBUG) 135 PRTDEN=( INDEX(KEYWRD,' DENS').NE.0 .AND. DEBUG) 136 PRTFOK=( INDEX(KEYWRD,'FOCK') .NE. 0 .AND. DEBUG) 137 PRTVEC=( INDEX(KEYWRD,'VECT') .NE. 0 .AND. DEBUG) 138 DEBUG=( INDEX(KEYWRD,'ITER') .NE. 0 ) 139C 140C INITIALIZE SOME LOGICALS AND CONSTANTS 141C 142 NEWDG =.FALSE. 143 PLCHEK=0.005D0 144 PL =1.D0 145 BSHIFT=0.D0 146 SHIFT=1.D0 147* 148* SCFCRT IS MACHINE-PRECISION DEPENDENT 149* 150 SCFCRT=1.D-4 151 ITRMAX = 200 152 NA2EL=NCLOSE 153 NA1EL=NALPHA+NOPEN 154 NB2EL=0 155 NB1EL=NBETA+NOPEN 156C 157C USE KEY-WORDS TO ASSIGN VARIOUS CONSTANTS 158C 159 IF(INDEX(KEYWRD,'FILL').NE.0) 160 1 IFILL=-READA(KEYWRD,INDEX(KEYWRD,'FILL')) 161 IF(INDEX(KEYWRD,'SHIFT').NE.0) 162 1 BSHIFT=-READA(KEYWRD,INDEX(KEYWRD,'SHIFT')) 163 IF(BSHIFT.NE.0)TEN=BSHIFT 164 IF(INDEX(KEYWRD,'ITRY').NE.0) 165 1 ITRMAX=READA(KEYWRD,INDEX(KEYWRD,'ITRY')) 166 CAMKIN=(INDEX(KEYWRD,'KING')+INDEX(KEYWRD,'CAMP') .NE. 0) 167 CI=(INDEX(KEYWRD,'MICROS')+INDEX(KEYWRD,'C.I.') .NE. 0) 168 OKPULY=.FALSE. 169 OKPULY=(INDEX(KEYWRD,'PULAY').NE.0) 170 UHF=(INDEX(KEYWRD,'UHF') .NE. 0) 171 SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0) 172 OKNEWD=ABS(BSHIFT) .LT. 0.001D0 173 IF(CAMKIN.AND.ABS(BSHIFT).GT.1.D-5) BSHIFT=4.44D0 174 TIMES=(INDEX(KEYWRD,'TIMES') .NE. 0) 175 TIMITR=(TIMES.AND.DEBUG) 176 FORCE=(INDEX(KEYWRD,'FORCE') .NE. 0) 177 ALLCON=(OKPULY.OR.CAMKIN) 178C 179C DO WE NEED A CAPPED ATOM CORRECTION? 180C 181 J=0 182 DO 10 I=1,NUMAT 183 10 IF(NAT(I).EQ.102)J=J+1 184 CAPPS=(J.GT.0) 185 IITER=1 186 TRANS=0.1D0 187 IF(INDEX(KEYWRD,'RESTART')+INDEX(KEYWRD,'OLDENS') 188 1 .NE. 0) THEN 189 IF(INDEX(KEYWRD,'OLDENS').NE.0) 190 1 OPEN(UNIT=10,FILE=GETNAM('FOR010'), 191 + STATUS='UNKNOWN',FORM='UNFORMATTED') 192 REWIND 10 193 READ(10)(PA(I),I=1,LINEAR) 194 IF( UHF) THEN 195 READ(10)(PB(I),I=1,LINEAR) 196 DO 20 I=1,LINEAR 197 POLD(I)=PA(I) 198 PBOLD(I)=PB(I) 199 20 P(I)=PA(I)+PB(I) 200 ELSE 201 DO 30 I=1,LINEAR 202 PB(I)=PA(I) 203 PBOLD(I)=PA(I) 204 POLD(I)=PA(I) 205 30 P(I)=PA(I)*2.D0 206 ENDIF 207 ELSE 208 NSCF=0 209 DO 40 I=1,LINEAR 210 P(I)=0.D0 211 PA(I)=0.D0 212 40 PB(I)=0.D0 213 W1=NA1EL/(NA1EL+1.D-6+NB1EL) 214 W2=1.D0-W1 215 IF(W1.LT.1.D-6)W1=0.5D0 216 IF(W2.LT.1.D-6)W2=0.5D0 217C 218C SLIGHTLY PERTURB THE DENSITY MATRIX IN CASE THE SYSTEM IS 219C TRAPPED IN A S**2 = 0 STATE. 220C 221 RANDOM=1.0D0 222 IF(UHF.AND.NA1EL.EQ.NB1EL) RANDOM=1.1D0 223 DO 50 I=1,NORBS 224 J=(I*(I+1))/2 225 P(J)=PDIAG(I) 226 PA(J)=P(J)*W1*RANDOM 227 RANDOM=1.D0/RANDOM 228 50 PB(J)=P(J)*W2*RANDOM 229 DO 60 I=1,LINEAR 230 PBOLD(I)=PB(I) 231 60 POLD(I)=PA(I) 232 ENDIF 233 HALFE=(NOPEN .NE. NCLOSE.AND.FRACT.NE.2.D0.AND.FRACT.NE.0.D0) 234C 235C DETERMINE THE SELF-CONSISTENCY CRITERION 236C 237 IF(INDEX(KEYWRD,'PREC') .NE. 0) 238 1 SCFCRT=SCFCRT*0.01D0 239 IF( INDEX(KEYWRD,'POLAR') + INDEX(KEYWRD,'NLLSQ') + 240 1 INDEX(KEYWRD,'SIGMA') .NE. 0) SCFCRT=SCFCRT*0.001D0 241 IF(FORCE) SCFCRT=SCFCRT*0.0001D0 242 IF(NOPEN-NCLOSE.GT.4) SCFCRT=SCFCRT*0.1D0 243 SCFCRT=MAX(SCFCRT,1.D-12) 244 IF(INDEX(KEYWRD,'POLAR').NE.0)SCFCRT=1.D-11 245C 246C THE USER CAN STATE THE SCF CRITERION, IF DESIRED. 247C 248 I=INDEX(KEYWRD,'SCFCRT') 249 IF(I.NE.0) THEN 250 SCFCRT=READA(KEYWRD,I) 251 WRITE(6,'('' SCF CRITERION ='',G14.4)')SCFCRT 252 IF(SCFCRT.LT.1.D-12) 253 1 WRITE(6,'(//2X,'' THERE IS A RISK OF INFINITE LOOPING WITH'', 254 2'' THE SCFCRT LESS THAN 1.D-12'')') 255 ELSE 256 IF(DEBUG)WRITE(6,'('' SCF CRITERION ='',G14.4)')SCFCRT 257 ENDIF 258 IF(.NOT.SCF1)LAST=0 259C 260C END OF INITIALIZATION SECTION. 261C 262 ELSEIF(FORCE.AND.NSCF.GT.0.AND..NOT.UHF)THEN 263C 264C RESET THE DENSITY MATRIX IF MECI HAS FORMED AN EXCITED STATE. THIS 265C PREVENTS THE SCF GETTING TRAPPED ON AN EXCITED STATE, PARTICULARLY 266C IF THE PULAY CONVERGER IS USED. 267C 268 DO 70 I=1,LINEAR 269 70 P(I)=2.D0*PA(I) 270 ENDIF 271C 272C INITIALIZATION OPERATIONS DONE EVERY TIME ITER IS CALLED 273C 274 MAKEA=.TRUE. 275 MAKEB=.TRUE. 276 IEMIN=0 277 IEMAX=0 278C 279C TURN OFF SHIFT IF NOT A FULL SCF. 280C 281 IF(.NOT.FULSCF) SHIFT=0.D0 282 IF(NEWDG) NEWDG=(ABS(BSHIFT).LT.0.001D0) 283 IF(LAST.EQ.1) NEWDG=.FALSE. 284C 285C SELF-CONSISTENCY CRITERIA: SELCON IS IN KCAL/MOL, PLTEST IS 286C A LESS IMPORTANT TEST TO MAKE SURE THAT THE SELCON TEST IS NOT 287C PASSED 'BY ACCIDENT' 288C IF GNORM IS LARGE, MAKE SELCON BIGGER 289C 290 SELCON=SCFCRT 291 IF(.NOT. FORCE .AND. .NOT. HALFE) THEN 292 IF(GNORM.GT.5.D0) SELCON=SCFCRT*GNORM*0.2D0 293 IF(GNORM.GT.200.D0) SELCON=SCFCRT*40.D0 294 ENDIF 295 PLTEST=0.05D0*SQRT(SELCON) 296C 297C SOMETIMES HEAT GOES SCF BUT DENSITY IS STILL FLUCTUATING IN UHF 298C IN WHICH CASE PAY LESS ATTENTION TO DENSITY MATRIX 299C 300 IF(NALPHA.NE.NBETA.AND.UHF)PLTEST=0.001D0 301 IF(DEBUG)THEN 302 WRITE(6,'('' SELCON, PLTEST'',3G16.7)')SELCON, PLTEST 303 ENDIF 304 IF(PRT1EL) THEN 305 WRITE(6,'(//10X,''ONE-ELECTRON MATRIX AT ENTRANCE TO ITER'')') 306 CALL VECPRT(H,NORBS) 307 ENDIF 308 IREDY=1 309 80 NITER=0 310 FRST=.TRUE. 311 IF(CAMKIN) THEN 312 MODEA=1 313 MODEB=1 314 ELSE 315 MODEA=0 316 MODEB=0 317 ENDIF 318 BFRST=.TRUE. 319********************************************************************** 320* * 321* * 322* START THE SCF LOOP HERE * 323* * 324* * 325********************************************************************** 326 INCITR=.TRUE. 327 90 INCITR=(MODEA.NE.3.AND.MODEB.NE.3) 328 IF(INCITR)NITER=NITER+1 329 IF(TIMITR)THEN 330 TITER=SECOND() 331 WRITE(6,*) 332 WRITE(6,'(A,F7.2)')' TIME FOR ITERATION:', TITER-TITER0 333 TITER0=TITER 334 ENDIF 335 IF(NITER.GT.ITRMAX-10.AND..NOT.ALLCON) THEN 336************************************************************************ 337* * 338* SWITCH ON ALL CONVERGERS * 339* * 340************************************************************************ 341 WRITE(6,'(//,'' ALL CONVERGERS ARE NOW FORCED ON'',/ 342 1 '' SHIFT=10, PULAY ON, CAMP-KING ON'',/ 343 2 '' AND ITERATION COUNTER RESET'',//)') 344 ALLCON=.TRUE. 345 BSHIFT=4.44D0 346 IREDY=-4 347 EOLD=100.D0 348 OKPULY=.TRUE. 349 NEWDG=.FALSE. 350 CAMKIN=(.NOT.HALFE) 351 GOTO 80 352 ENDIF 353************************************************************************ 354* * 355* MAKE THE ALPHA FOCK MATRIX * 356* * 357************************************************************************ 358 IF(ABS(SHIFT).GT.1.D-10.AND.BSHIFT .NE. 0.D0) THEN 359 L=0 360 IF(NITER.GT.1)THEN 361 IF(NEWDG.AND..NOT.(HALFE.OR.CAMKIN))THEN 362C 363C SHIFT WILL APPLY TO THE VIRTUAL ENERGY LEVELS USED IN THE 364C PSEUDODIAGONALIIZATION. IF DIFF IS -VE, GOOD, THEN LOWER THE 365C HOMO-LUMO GAP BY 0.1EV, OTHERWISE INCREASE IT. 366 IF(DIFF.GT.0)THEN 367 SHIFT=1.D0 368C 369C IF THE PSEUDODIAGONALIZATION APPROXIMATION -- THAT THE WAVEFUNCTION 370C IS ALMOST STABLE -- IS INVALID, TURN OFF NEWDG 371 IF(DIFF.GT.1)NEWDG=.FALSE. 372 ELSE 373 SHIFT=-0.1D0 374 ENDIF 375 ELSE 376 SHIFT=TEN+EIGS(IHOMO+1)-EIGS(IHOMO)+SHIFT 377 ENDIF 378 IF(DIFF.GT.0.D0) THEN 379 IF(SHIFT.GT.4.D0)SHFMAX=4.5D0 380 IF(SHIFT.GT.SHFMAX)SHFMAX=MAX(SHFMAX-0.5D0,0.D0) 381 ENDIF 382C 383C IF SYSTEM GOES UNSTABLE, LIMIT SHIFT TO THE RANGE -INFINITY - SHFMAX 384C BUT IF SYSTEM IS STABLE, LIMIT SHIFT TO THE RANGE -INFINITY - +20 385C 386 SHIFT=MAX(-20.D0,MIN(SHFMAX,SHIFT)) 387 IF(ABS(SHIFT-SHFMAX).LT.1.D-5)SHFMAX=SHFMAX+0.01D0 388C 389C THE CAMP-KING AND PULAY CONVERGES NEED A CONSTANT SHIFT. 390C IF THE SHIFT IS ALLOWED TO VARY, THESE CONVERGERS WILL NOT 391C WORK PROPERLY. 392C 393 IF(OKPULY.OR.ABS(BSHIFT-4.44D0).LT.1.D-5)THEN 394 SHIFT=-8.D0 395 IF(NEWDG) SHIFT=0 396 ENDIF 397 IF(UHF)THEN 398 IF(NEWDG.AND..NOT.(HALFE.OR.CAMKIN))THEN 399 SHIFTB=TEN-TENOLD 400 ELSE 401 SHIFTB=TEN+EIGS(IHOMOB+1)-EIGS(IHOMOB)+SHIFTB 402 ENDIF 403 IF(DIFF.GT.0.D0)SHIFTB=MIN(4.D0,SHIFTB) 404 SHIFTB=MAX(-20.D0,MIN(SHFMAX,SHIFTB)) 405 IF(OKPULY.OR.ABS(BSHIFT-4.44D0).LT.1.D-5)THEN 406 SHIFTB=-8.D0 407 IF(NEWDG)SHIFT=0 408 ENDIF 409 DO 100 I=IHOMOB+1,NORBS 410 100 EIGB(I)=EIGB(I)+SHIFTB 411 ELSE 412 ENDIF 413 ENDIF 414 TENOLD=TEN 415 IF(PL.GT.PLCHEK)THEN 416 SHFTBO=SHIFTB 417 SHFTO=SHIFT 418 ELSE 419 SHIFTB=SHFTBO 420 SHIFT=SHFTO 421 ENDIF 422 DO 110 I=IHOMO+1,NORBS 423 110 EIGS(I)=EIGS(I)+SHIFT 424 DO 130 I=1,NORBS 425 DO 120 J=1,I 426 L=L+1 427 120 F(L)=H(L)+SHIFT*PA(L) 428 130 F(L)=F(L)-SHIFT 429 ELSEIF (I.EQ.77.AND.LAST.EQ.0.AND.NITER.LT.2.AND.FULSCF)THEN 430C 431C SLIGHTLY PERTURB THE FOCK MATRIX IN CASE THE SYSTEM IS 432C TRAPPED IN A METASTABLE EXCITED ELECTRONIC STATE 433C 434 RANDOM=0.001D0 435 DO 140 I=1,LINEAR 436 RANDOM=-RANDOM 437 140 F(I)=H(I)+RANDOM 438 ELSE 439 DO 150 I=1,LINEAR 440 150 F(I)=H(I) 441 ENDIF 442 160 CONTINUE 443 IF(TIMITR)THEN 444 T0=SECOND() 445 WRITE(6,'(A,F7.2)')' LOAD FOCK MAT. INTEGRAL',T0-TITER0 446 ENDIF 447C# CALL TIMER('BEFORE FOCK2') 448 CALL FOCK2(F,P,PA,W, WJ, WK,NUMAT,NAT,NFIRST,NMIDLE,NLAST) 449C# CALL TIMER('AFTER FOCK2') 450C# CALL TIMER('BEFORE FOCK1') 451 CALL FOCK1(F,P,PA,PB) 452C# CALL TIMER('AFTER FOCK1') 453 IF(TIMITR)THEN 454 T0=SECOND() 455 TF1=TF1+T0-T1 456 WRITE(6,'(2(A,F7.2))')' FOCK1:',T0-T1,'INTEGRAL:',T0-TITER0 457 ENDIF 458************************************************************************ 459* * 460* MAKE THE BETA FOCK MATRIX * 461* * 462************************************************************************ 463 IF (UHF) THEN 464 IF(SHIFTB .NE. 0.D0) THEN 465 L=0 466 DO 180 I=1,NORBS 467 DO 170 J=1,I 468 L=L+1 469 170 FB(L)=H(L)+SHIFTB*PB(L) 470 180 FB(L)=FB(L)-SHIFTB 471 ELSEIF (RAND.AND.LAST.EQ.0.AND.NITER.LT.2.AND.FULSCF)THEN 472 RANDOM=0.001D0 473 DO 190 I=1,LINEAR 474 RANDOM=-RANDOM 475 190 FB(I)=H(I)+RANDOM 476 ELSE 477 DO 200 I=1,LINEAR 478 200 FB(I)=H(I) 479 ENDIF 480 CALL FOCK2(FB,P,PB,W, WJ, WK,NUMAT,NAT,NFIRST,NMIDLE,NLAST) 481 CALL FOCK1(FB,P,PB,PA) 482 ENDIF 483 IF( .NOT. FULSCF) GOTO 380 484 IF(PRTFOK) THEN 485 WRITE(6,210)NITER 486 210 FORMAT(' FOCK MATRIX ON ITERATION',I3) 487 CALL VECPRT (F,NORBS) 488 ENDIF 489C 490C CODE THE FOLLOWING LINE IN PROPERLY SOMETIME 491C THIS OPERATION IS BELIEVED TO GIVE RISE TO A BETTER FOCK MATRIX 492C THAN THE CONVENTIONAL GUESS. 493C 494 IF(IRRR.EQ.0)THEN 495 DO 220 I=1,NORBS 496 220 F((I*(I+1))/2)=F((I*(I+1))/2)*0.5D0 497 IRRR=2 498 ENDIF 499************************************************************************ 500* * 501* CALCULATE THE ENERGY IN KCAL/MOLE * 502* * 503************************************************************************ 504 IF (NITER .GE. ITRMAX) THEN 505 IF(DIFF.LT.1.D-3.AND.PL.LT.1.D-4.AND..NOT.FORCE)THEN 506 WRITE(6,'('' """""""""""""""UNABLE TO ACHIEVE SELF-CONSISTEN 507 1CE, JOB CONTINUING'')') 508 GOTO 380 509 ENDIF 510 IF(MINPRT)WRITE (6,230) 511 230 FORMAT (//10X,'"""""""""""""UNABLE TO ACHIEVE SELF-CONSISTENCE' 512 1,/) 513 WRITE (6,240) DIFF,PL 514 240 FORMAT (//,10X,'DELTAE= ',E12.4,5X,'DELTAP= ',E12.4,///) 515C *** here we failed to calculate a valid energy, but we don't want to close the whole program either. 516C *** instead of calling STOP, continue like in the above case where GOTO 380 is called... 517 GOTO 380 518C IFLEPO=9 519C IITER=2 520C CALL WRITMO(TIME0,ESCF) 521C STOP 522 ENDIF 523 EE=HELECT(NORBS,PA,H,F) 524 IF(UHF)THEN 525 EE=EE+HELECT(NORBS,PB,H,FB) 526 ELSE 527 EE=EE*2.D0 528 ENDIF 529 IF(CAPPS)EE=EE+CAPCOR(NAT,NFIRST,NLAST,NUMAT,P,H) 530 IF(BSHIFT.NE.0.D0) 531 1SCORR=SHIFT*(NOPEN-NCLOSE)*23.061D0*0.25D0*(FRACT*(2.D0-FRACT)) 532 ESCF=(EE+ENUCLR)*23.061D0+ATHEAT+SCORR 533 IF(INCITR)THEN 534 DIFF=ESCF-EOLD 535 IF(DIFF.GT.0)THEN 536 TEN=TEN-1.D0 537 ELSE 538 TEN=TEN*0.975D0+0.05D0 539 ENDIF 540C 541C MAKE SURE SELF-CONSISTENCY TEST IS NOT MORE STRINGENT THAN THE 542C COMPUTER CAN HANDLE 543C 544 SELLIM=MAX(SELCON,EPS*MAX(ABS(EE),1.D0)) 545C 546C SCF TEST: CHANGE IN HEAT OF FORMATION IN KCAL/MOL SHOULD BE 547C LESS THAN SELLIM. THE OTHER TESTS ARE SAFETY MEASURES 548C 549 IF(.NOT.(NITER.GT.4.AND.(PL.EQ.0.D0.OR.PL.LT.PLTEST.AND. 550 1 ABS(DIFF).LT.SELLIM) .AND. READY)) GOTO 270 551************************************************************************ 552* * 553* SELF-CONSISTENCY TEST, EXIT MODE FROM ITERATIONS * 554* * 555************************************************************************ 556 250 IF (ABS(SHIFT) .LT. 1.D-10) GOTO 380 557 SHIFT=0.D0 558 SHIFTB=0.D0 559 DO 260 I=1,LINEAR 560 260 F(I)=H(I) 561 MAKEA=.TRUE. 562 MAKEB=.TRUE. 563 GOTO 160 564 270 CONTINUE 565*********************************************************************** 566*********************************************************************** 567 IF(LIMSCF.AND.EMIN.NE.0.D0.AND..NOT.(CI.OR.HALFE))THEN 568C 569C THE FOLLOWING TESTS ARE INTENDED TO ALLOW A FAST EXIT FROM ITER 570C IF THE RESULT IS 'GOOD ENOUGH' FOR THE CURRENT STEP IN THE GEOMETRY 571C OPTIMIZATION 572C 573 IF(ESCF.LT.EMIN)THEN 574C 575C THE ENERGY IS LOWER THAN THE PREVIOUS MINIMUM. NOW CHECK THAT 576C IT IS CONSISTENTLY LOWER. 577C 578 IEMAX=0 579 IEMIN=MIN(5,IEMIN+1) 580 DO 280 I=2,IEMIN 581 280 ESCF0(I-1)=ESCF0(I) 582 ESCF0(IEMIN)=ESCF 583C 584C IS THE DIFFERENCE IN ENERGY BETWEEN TWO ITERATIONS LESS THAN 5% 585C OF THE ENERGY GAIN FOR THIS GEOMETRY RELATIVE TO THE PREVIOUS 586C MINIMUM. 587C 588 IF(IEMIN.GT.3)THEN 589 DO 290 I=2,IEMIN 590 IF(ABS(ESCF0(I)-ESCF0(I-1)).GT.0.05D0*(EMIN-ESCF)) 591 1GOTO 320 592 290 CONTINUE 593C 594C IS GOOD ENOUGH -- RAPID EXIT 595C 596 IF(DEBUG) WRITE(6,*) 597 1' RAPID EXIT BECAUSE ENERGY IS CONSISTENTLY LOWER' 598 GOTO 250 599 ENDIF 600 ELSE 601C 602C THE ENERGY HAS RISEN ABOVE THAT OF THE PREVIOUS MINIMUM. 603C WE NEED TO CHECK WHETHER THIS IS A FLUKE OR IS THIS REALLY 604C A BAD GEOMETRY. 605C 606 IEMIN=0 607 IEMAX=MIN(5,IEMAX+1) 608 DO 300 I=2,IEMAX 609 300 ESCF0(I-1)=ESCF0(I) 610 ESCF0(IEMAX)=ESCF 611C 612C IS THE DIFFERENCE IN ENERGY BETWEEN TWO ITERATIONS LESS THAN 5% 613C OF THE ENERGY LOST FOR THIS GEOMETRY RELATIVE TO THE PREVIOUS 614C MINIMUM. 615C 616 IF(IEMAX.GT.3)THEN 617 DO 310 I=2,IEMAX 618 IF(ABS(ESCF0(I)-ESCF0(I-1)).GT.0.05D0*(ESCF-EMIN)) 619 1GOTO 320 620 310 CONTINUE 621C 622C IS GOOD ENOUGH -- RAPID EXIT 623C 624 IF(DEBUG) WRITE(6,*) 625 1' RAPID EXIT BECAUSE ENERGY IS CONSISTENTLY HIGHER' 626 GOTO 250 627 ENDIF 628 ENDIF 629 ENDIF 630 320 READY=(IREDY.GT.0.AND.(ABS(DIFF).LT.SELLIM*10.D0.OR.PL.EQ.0.D0) 631 1) 632 IREDY=IREDY+1 633 ENDIF 634 IF(PRTPL.OR.DEBUG.AND.NITER.GT.ITRMAX-20) THEN 635 IF(ABS(ESCF).GT.99999.D0) ESCF=SIGN(9999.D0,ESCF) 636 IF(ABS(DIFF).GT.9999.D0)DIFF=0.D0 637 IF(INCITR) 638 1 WRITE(6,'('' ITERATION'',I3,'' PLS='',2E10.3,'' ENERGY '', 639 2F14.7,'' DELTAE'',F13.7)')NITER,PL,PLB,ESCF,DIFF 640 close (6) 641C ***** Modified by Jiro Toyoda at 1994-05-25 ***** 642C OPEN(UNIT=6,FILE=GETNAM('FOR006'),ACCESS='APPEND') 643C *** exactly why do we want to open unit 6??? it's already open?!?!?! 644C *** we also remove this because we want use STDOUT for output... 645C OPEN(UNIT=6,FILE=GETNAM('FOR006')) 646C 9990 read (6,'()',end=9999) 647C goto 9990 648C 9999 continue 649C ***************************** at 1994-05-25 ***** 650 ENDIF 651 IF(INCITR)EOLD=ESCF 652************************************************************************ 653* * 654* INVOKE THE CAMP-KING CONVERGER * 655* * 656************************************************************************ 657 IF(NITER.GT.2 .AND. CAMKIN .AND. MAKEA) 658 1CALL INTERP(NORBS,NA1EL,NORBS-NA1EL, MODEA, ESCF/23.061D0, 659 2F, C, AR1, AR2, AR3, AR4, AR1) 660 MAKEB=.FALSE. 661 IF(MODEA.EQ.3)GOTO 340 662 MAKEB=.TRUE. 663 IF(TIMITR)THEN 664 T0=SECOND() 665 WRITE(6,'(2(A,F7.2))')' ADJUST DAMPER INTEGRAL',T0-TITER0 666 ENDIF 667 IF( NEWDG ) THEN 668************************************************************************ 669* * 670* INVOKE PULAY'S CONVERGER * 671* * 672************************************************************************ 673 IF(OKPULY.AND.MAKEA.AND.IREDY.GT.1) 674 1CALL PULAY(F,PA,NORBS,POLD,POLD2,POLD3,JALP,IALP,MPACK,FRST,PL) 675************************************************************************ 676* * 677* DIAGONALIZE THE ALPHA OR RHF SECULAR DETERMINANT * 678* WHERE POSSIBLE, USE THE PULAY-STEWART METHOD, OTHERWISE USE BEPPU'S * 679* * 680************************************************************************ 681 IF (HALFE.OR.CAMKIN) THEN 682 CALL HQRII(F,NORBS,NORBS,EIGS,C) 683 ELSE 684C# CALL TIMER('BEFORE DIAG') 685 CALL DIAG (F,C,NA1EL,EIGS,NORBS,NORBS) 686C# CALL TIMER('AFTER DIAG') 687 ENDIF 688 ELSE 689C# CALL TIMER('BEFORE HQRII') 690 CALL HQRII(F,NORBS,NORBS,EIGS,C) 691C# CALL TIMER('AFTER HQRII') 692 IF(TIMITR)THEN 693 T1=SECOND() 694 WRITE(6,'(2(A,F7.2))')' HQRII:',T1-T0,' INTEGRAL',T1-TITER0 695 ENDIF 696 ENDIF 697 J=1 698 IF(PRTVEC) THEN 699 J=1 700 IF(UHF)J=2 701 WRITE(6,'(//10X,A, 702 1'' EIGENVECTORS AND EIGENVALUES ON ITERATION'',I3)') 703 2 ABPRT(J),NITER 704 CALL MATOUT(C,EIGS,NORBS,NORBS,NORBS) 705 ELSE 706 IF (PRTEIG) WRITE(6,330)ABPRT(J),NITER,(EIGS(I),I=1,NORBS) 707 ENDIF 708 330 FORMAT(10X,A,' EIGENVALUES ON ITERATION',I3,/10(6G13.6,/)) 709 340 IF(IFILL.NE.0)CALL SWAP(C,NORBS,NORBS,NA2EL,IFILL) 710************************************************************************ 711* * 712* CALCULATE THE ALPHA OR RHF DENSITY MATRIX * 713* * 714************************************************************************ 715 IF(UHF)THEN 716 CALL DENSIT( C,NORBS, NORBS, NA2EL,NA1EL, FRACT, PA, 1) 717 IF(MODEA.NE.3.AND..NOT. (NEWDG.AND.OKPULY)) 718 1 CALL CNVG(PA, POLD, POLD2, NORBS, NITER, PL) 719 ELSE 720C# CALL TIMER('BEFORE DENSIT') 721 CALL DENSIT( C,NORBS, NORBS, NA2EL,NA1EL, FRACT, P, 1) 722C# CALL TIMER('AFTER DENSIT') 723 IF(MODEA.NE.3.AND..NOT. (NEWDG.AND.OKPULY))THEN 724C# CALL TIMER('BEFORE CNVG') 725 CALL CNVG(P, POLD, POLD2, NORBS, NITER, PL) 726C# CALL TIMER('AFTER CNVG') 727 ENDIF 728 ENDIF 729************************************************************************ 730* * 731* UHF-SPECIFIC CODE * 732* * 733************************************************************************ 734 IF( UHF )THEN 735************************************************************************ 736* * 737* INVOKE THE CAMP-KING CONVERGER * 738* * 739************************************************************************ 740 IF(NITER.GT.2 .AND. CAMKIN .AND. MAKEB ) 741 1CALL INTERP(NORBS,NB1EL,NORBS-NB1EL, MODEB, ESCF/23.061D0, 742 2FB, CBETA, BR1, BR2, BR3, BR4, BR1) 743 MAKEA=.FALSE. 744 IF(MODEB.EQ.3) GOTO 350 745 MAKEA=.TRUE. 746 IF( NEWDG ) THEN 747************************************************************************ 748* * 749* INVOKE PULAY'S CONVERGER * 750* * 751************************************************************************ 752 IF( OKPULY.AND.MAKEB.AND.IREDY.GT.1) 753 1CALL PULAY(FB,PB,NORBS,PBOLD,PBOLD2, 754 2PBOLD3,JBET,IBET,MPACK,BFRST,PLB) 755************************************************************************ 756* * 757* DIAGONALIZE THE ALPHA OR RHF SECULAR DETERMINANT * 758* WHERE POSSIBLE, USE THE PULAY-STEWART METHOD, OTHERWISE USE BEPPU'S * 759* * 760************************************************************************ 761 IF (HALFE.OR.CAMKIN) THEN 762 CALL HQRII(FB,NORBS,NORBS,EIGB,CBETA) 763 ELSE 764 CALL DIAG (FB,CBETA,NB1EL,EIGB,NORBS,NORBS) 765 ENDIF 766 ELSE 767 CALL HQRII(FB,NORBS,NORBS,EIGB,CBETA) 768 ENDIF 769 IF(PRTVEC) THEN 770 WRITE(6,'(//10X,A,'' EIGENVECTORS AND EIGENVALUES ON '', 771 1''ITERATION'',I3)')ABPRT(3),NITER 772 CALL MATOUT(CBETA,EIGB,NORBS,NORBS,NORBS) 773 ELSE 774 IF (PRTEIG) WRITE(6,330)ABPRT(3),NITER,(EIGB(I),I=1,NORBS) 775 ENDIF 776************************************************************************ 777* * 778* CALCULATE THE BETA DENSITY MATRIX * 779* * 780************************************************************************ 781 350 CALL DENSIT( CBETA,NORBS, NORBS, NB2EL, NB1EL, FRACT, PB, 1) 782 IF( .NOT. (NEWDG.AND.OKPULY)) 783 1CALL CNVG(PB, PBOLD, PBOLD2, NORBS, NITER, PLB) 784 ENDIF 785************************************************************************ 786* * 787* CALCULATE THE TOTAL DENSITY MATRIX * 788* * 789************************************************************************ 790 IF(UHF) THEN 791 DO 360 I=1,LINEAR 792 360 P(I)=PA(I)+PB(I) 793 ELSE 794 DO 370 I=1,LINEAR 795 PA(I)=P(I)*0.5D0 796 370 PB(I)=PA(I) 797 ENDIF 798 IF(PRTDEN) THEN 799 WRITE(6,'('' DENSITY MATRIX ON ITERATION'',I4)')NITER 800 CALL VECPRT (P,NORBS) 801 ENDIF 802 OKNEWD=(PL.LT.SELLIM .OR. OKNEWD) 803 NEWDG=(PL.LT.TRANS .AND. OKNEWD .OR. NEWDG) 804 IF(PL.LT.TRANS*0.3333D0)OKNEWD=.TRUE. 805 GO TO 90 806********************************************************************** 807* * 808* * 809* END THE SCF LOOP HERE * 810* NOW CALCULATE THE ELECTRONIC ENERGY * 811* * 812* * 813********************************************************************** 814* SELF-CONSISTENCE ACHEIVED. 815* 816 380 EE=HELECT(NORBS,PA,H,F) 817 IF(UHF) THEN 818 EE=EE+HELECT(NORBS,PB,H,FB) 819 ELSE 820 EE=EE*2.D0 + 821 1SHIFT*(NOPEN-NCLOSE)*23.061D0*0.25D0*(FRACT*(2.D0-FRACT)) 822 ENDIF 823 IF(CAPPS)EE=EE+CAPCOR(NAT,NFIRST,NLAST,NUMAT,P,H) 824C 825C NORMALLY THE EIGENVALUES ARE INCORRECT BECAUSE THE 826C PSEUDODIAGONALIZATION HAS BEEN USED. IF THIS 827C IS THE LAST SCF, THEN DO AN EXACT DIAGONALIZATION 828 IF( NSCF.EQ.0 .OR. LAST.EQ.1 .OR. CI .OR. HALFE ) THEN 829C 830C PUT F AND FB INTO POLD IN ORDER TO NOT DESTROY F AND FB 831C AND DO EXACT DIAGONALISATIONS 832 DO 390 I=1,LINEAR 833 390 POLD(I)=F(I) 834 CALL HQRII(POLD,NORBS,NORBS,EIGS,C) 835 IF(UHF) THEN 836 DO 400 I=1,LINEAR 837 400 POLD(I)=FB(I) 838 CALL HQRII(POLD,NORBS,NORBS,EIGB,CBETA) 839 DO 410 I=1,LINEAR 840 410 POLD(I)=PA(I) 841 ELSE 842 DO 420 I=1,LINEAR 843 420 POLD(I)=P(I) 844 ENDIF 845 IF(CI.OR.HALFE) THEN 846C# CALL TIMER('BEFORE MECI') 847 SUM=MECI(EIGS,C) 848C# CALL TIMER('AFTER MECI') 849 EE=EE+SUM 850 IF(PRTPL)THEN 851 ESCF=(EE+ENUCLR)*23.061D0+ATHEAT 852 WRITE(6,'(27X,''AFTER MECI, ENERGY '',F14.7)')ESCF 853 ENDIF 854 ENDIF 855 ENDIF 856 NSCF=NSCF+1 857 IF(DEBUG)WRITE(6,'('' NO. OF ITERATIONS ='',I6)')NITER 858C IF(FORCE) SCFCRT=1.D-5 859 IF(ALLCON.AND.ABS(BSHIFT-4.44D0).LT.1.D-7)THEN 860 CAMKIN=.FALSE. 861 ALLCON=.FALSE. 862 NEWDG=.FALSE. 863 BSHIFT=-10.D0 864 OKPULY=.FALSE. 865 ENDIF 866 SHIFT=1.D0 867 IF(EMIN.EQ.0.D0)THEN 868 EMIN=ESCF 869 ELSE 870 EMIN=MIN(EMIN,ESCF) 871 ENDIF 872 RETURN 873 END 874