program umat implicit none integer*4 ni,nc,i,j,ia,ix,ic,nz,IERR,m,ij,il,jj,N1,N2,ii,ne,jf,kk, 1nat,nb real*8 t,am,tol,oi,nii,spii,t2 real*8,allocatable::B(:,:),G(:,:),K(:,:),L(:),NK(:,:) character*160 s160 tol=1.0d-2 t2=1.0d-4 write(6,601) 601 format(' Non-redundant combination of interanl coordinates',/, 1 ' Input: OLD.FILE.UMA',/, 1 ' OLD.B.MAT',/, 1 ' Output: FILE.UMA') open(9,file='OLD.B.MAT',status='old') read(9,*)ni read(9,*)nc read(9,*) allocate(B(ni,nc)) B=0.0d0 do 1 i=1,ni read(9,*) read(9,*)m do 1 j=1,m 1 read(9,*)ia,ix,B(i,ix+3*(ia-1)) close(9) write(6,605)ni,nc 605 format(i5,' internal,',i5,' Cartesian') allocate(G(ni,ni),K(ni,ni),L(ni)) G=0.0d0 do 2 i=1,ni do 2 j=1,ni do 2 ic=1,nc 2 G(i,j)=G(i,j)+B(i,ic)*B(j,ic) call TRED12(ni,G,K,L,2,IERR) nz=0 do 3 i=1,ni 3 if(dabs(L(i)).gt.tol)nz=nz+1 write(6,*)nz,' non-zero:' if(nz.lt.20)write(6,602)(L(i),i=1,nz) 602 format(5(10f10.3)) write(6,*) if(nz.lt.20)call wk(K,nz,ni) write(6,*)' Conserve N1 N2:' read(5,*)N1,N2 c put N1..N2 at the beginning do 52 i=N1,N2 ij=0 am=0.0d0 do 51 j=1+i-N1,nz if(dabs(K(i,j)).gt.am)then am=dabs(K(i,j)) ij=j endif 51 continue if(ij.eq.0.or.am.lt.tol)then write(6,*)'no solution' stop endif c N1 to line 1, N1+2 to line 2 etc il=i-N1+1 do 52 j=1,ni t=K(j,ij) K(j,ij)=K(j,il) 52 K(j,il)=t write(6,*)'reordered' c zero-out coupling terms: jf=0 allocate(NK(ni,ni)) NK=K do 5 j=1,ni do 55 i=1,N2-N1+1 if((j.lt.N1.or.j.gt.N2).and.dabs(K(j,i)).gt.t2)then if(jf.eq.0)jf=j ij=0 am=0.0d0 ne=nz if(jf.ne.j)ne=N2-N1+1 do 53 ii=1,ne if(dabs(K(j,ii)).gt.am.and.dabs(K(j,ii)).gt.tol.and.i.ne.ii)then c check that lines i and ii are not linearly dependent: oi=0.0d0 nii=0.0d0 spii=0.0d0 do 57 kk=1,ni oi =oi +K(kk,i )*K(kk,i) nii =nii +K(kk,ii)*K(kk,ii) 57 spii=spii+K(kk,i )*K(kk,ii) if(dabs(dsqrt(oi*nii)-dabs(spii)).gt.tol)then am=dabs(K(j,ii)) ij=ii endif endif 53 continue if(ij.eq.0)then write(6,609)i,j,K(j,i),ne 609 format(' i j:',2i3,' K:',f7.3, ' ne:',i3) write(6,*) call wk(K,nz,ni) write(6,*)'impossible to get rid of the coupling' stop endif do 54 jj=1,ni 54 NK(jj,i)=NK(jj,i)-K(jj,ij)*K(j,i)/K(j,ij) endif 55 continue K=NK 5 continue c delete mixed terms do 58 i=N2-N1+2,nz do 58 j=N1,N2 58 K(j,i)=0.0d0 do 59 i=1,N2-N1+1 do 591 j=1,N2-N1+1 591 K(j,i)=0.0d0 59 K(i,i)=1.0d0 write(6,*)'deleted coupling:' call wk(K,nz,ni) open(8,file='OLD.FILE.UMA',status='old') open(9,file='FILE.UMA') read(8,800)s160 800 format(a160) call wr(s160) read(8,800)s160 call wr(s160) read(8,800)s160 call wr(s160) read(s160,*)nat,nb do 60 i=1,nat read(8,800)s160 60 call wr(s160) do 61 i=1,nb read(8,800)s160 call wr(s160) read(s160,*)m if(m.eq.4)then read(8,800)s160 call wr(s160) read(8,800)s160 call wr(s160) endif 61 continue close(8) write(9,*)nz do 100 i=1,nz m=0 do 101 j=1,ni 101 if(dabs(K(j,i)).gt.tol)m=m+1 write(6,607)m,i write(9,607)m,i 607 format(2i5) do 100 j=1,ni if(dabs(K(j,i)).gt.tol)then write(6,608)j,K(j,i) write(9,608)j,K(j,i) 608 format(5x,i4,f7.3) endif 100 continue close(9) 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 subroutine wk(K,nz,ni) implicit none integer*4 nz,ni,i,j real*8 K(ni,ni) do 4 i=1,nz 4 write(6,603)(K(j,i),j=1,ni) 603 format(20f7.3) return end subroutine wr(s160) implicit none integer*4 ne,i character*160 s160 do 60 ne=len(s160),1,-1 60 if(s160(ne:ne).ne.' ')goto 61 61 do 62 i=1,ne 62 write(9,801)s160(i:i) 801 format(a1,$) write(9,*) return end