program checkfock implicit none integer*4 ln,n,io,iu,ia,i,j,k,IERR real*8 tol real*8,allocatable::c(:,:),F(:,:),s(:,:),E(:), 1si(:,:),U(:,:),L(:),fu(:,:),ufu(:,:),luful(:,:), 1D(:,:),EN(:),CN(:,:) character*80 s80 ln=0 n=0 io=0 iu=0 ia=0 open(9,file='G.OUT') 1 read(9,900,end=99,err=99)s80 900 format(a80) ln=ln+1 if(s80(8:22).eq.'basis functions')read(s80(1:7),*)n if(s80(2:16).eq.'*** Overlap ***'.and.n.ne.0)then write(6,900)s80 io=io+1 if(io.eq.1)allocate(s(n,n),si(n,n)) call rm(9,S,n,0) tol=1.0d-10 call INV(s,sI,N,TOL) endif if(s80(2:21).eq.'Fock matrix (alpha):'.and.n.ne.0)then write(6,900)s80 iu=iu+1 if(iu.eq.1)allocate(f(n,n)) call rm(9,F,n,0) endif if(s80(2:11).eq.'Alpha MOs:'.and.n.ne.0)then write(6,900)s80 ia=ia+1 if(ia.eq.1)allocate(c(n,n),e(n),L(n),U(n,N),fu(n,n), 1 ufu(n,n),luful(n,n),D(n,n),en(n),cn(n,n)) call ro(9,C,n,e) IERR=0 call TRED12(N,s,U,L,2,IERR) write(6,*)'s-EVs:' write(6,600)(L(i),i=1,n) 600 format(5f10.5) do 2 i=1,n do 2 J=1,n fu(i,j)=0.0d0 do 2 k=1,n 2 fu(i,j)=fu(i,j)+f(i,k)*u(k,J) do 3 I=1,n do 3 J=1,n ufu(i,j)=0.0d0 do 3 k=1,n 3 ufu(i,j)=ufu(i,j)+u(k,I)*fu(k,J) do 4 I=1,n do 4 J=1,n 4 luful(i,j)=ufu(i,j)/dsqrt(L(i)*L(j)) call TRED12(N,luful,D,EN,2,IERR) write(6,*)'new Es:' write(6,600)(EN(i),i=1,n) do 5 i=1,n do 5 J=1,n cn(i,J)=0.0d0 do 5 k=1,n 5 cn(i,J)=cn(i,J)+D(k,J)/dsqrt(L(k))*U(i,k) call wmtr(6,cn,n,'New Cs:',1) endif goto 1 99 close(9) write(6,*)ln,' lines' end subroutine wmtr(io,A,n,st,ic) c c ic=0 .. writes triangle of a supposedly symmetric matrix c ic=1 .. writes it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(n,n) character*(*) st write(io,*)st write(io,*) N1=1 1 N3=MIN(N1+4,N) WRITE(io,18)(I,I=N1,N3) 18 FORMAT(4X,5I14) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=n3 130 WRITE(io,17)LN,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) write(io,*) return end subroutine ro(io,A,n,e) c reads energies and orbitals IMPLICIT none INTEGER*4 n,N1,N3,ln,io,j real*8 A(n,n),e(*) N1=1 1 N3=MIN(N1+4,N) read(io,*) read(io,100)(e(J),J=N1,N3) 100 format(21x,5f10.5) DO 130 LN=1,N 130 read(io,100)(A(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.N)GOTO 1 return end subroutine rm(io,A,n,ic) c c ic=0 .. reads triangle of a supposedly symmetric matrix c ic=1 .. read it all c IMPLICIT none INTEGER*4 n,ic,N1,N3,lnst,ln,iend,i,io,j real*8 A(n,n) N1=1 1 N3=MIN(N1+4,N) read(io,*) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=n3 130 read(io,*)A(LN,N1),(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 if(ic.eq.0)then do 2 i=1,N do 2 j=1,i 2 A(J,i)=A(i,J) endif return end SUBROUTINE INV(A,AI,N,TOL) IMPLICIT none integer*4 N,oo,ii,jj,iw,io,kk,i2,j2 REAL*8 TOL, A(N,N),AI(N,N),w real*8, allocatable ::E(:,:) allocate(E(N,2*N)) DO 1 ii=1,N DO 1 jj=1,N e(ii,jj)=a(ii,jj) E(II,JJ+N)=0.0D0 1 if (ii.EQ.jj)e(ii,jj+N)=1.0D0 DO 2 ii=1,N-1 iw=ii if (DABS(e(iw,iw)).LT.TOL) then DO 3 io=iw+1,N oo=IO if (io.GT.N)oo=io-N 3 if (DABS(e(oo,iw)).GE.TOL) goto 11 write(6,*)'Inverse cannot be done' stop 11 CONTINUE DO 4 kk=1, 2*N w=e(iw,kk) e(iw,kk)=e(oo,kk) 4 e(oo,kk)=w ENDIF DO 2 jj=ii+1,N DO 6 kk=ii+1, 2*N 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii) 2 e(jj,ii)=0.0D0 DO 7 ii=1, N-1 i2=N-ii+1 DO 7 jj=1,i2-1 j2=i2-jj DO 9 kk=1, N 9 e(j2,kk+N)=e(j2,kk+N)-e(i2,kk+N)*e(j2,i2)/e(i2,i2) 7 e(j2,i2)=0.0D0 DO 10 ii=1,N DO 10 jj=1,N 10 AI(ii,jj)=e(ii,jj+N)/e(ii,ii) deallocate(E) RETURN END SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END