PROGRAM CSCALE IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*300,MX=MX3/3) CHARACTER*80 f1,f2,f3,f4 DIMENSION FCAR1(MX3,MX3),FCAR2(MX3,MX3), 1FCAR3(MX3,MX3),FCAR4(MX3,MX3), 2b1(MX3,MX3),b2(MX3,MX3),d1(MX3),d2(MX3),sf(MX3) C WRITE(6,60000) 60000 FORMAT(' Force field perturbational scaling',/,/, 1 ' Unperturbed ----> Perturbed',/, 2 ' Better Un ----> Better Per',/) write(6,*)'Unperturbed FF file name:' read(5,'(a)')f1 write(6,*)'Number of atoms:' read(5,*)nat1 write(6,*) na1=nat1*3 write(6,*)'Perturbed FF file name:' read(5,'(a)')f2 write(6,*)'Number of atoms:' read(5,*)nat2 na2=nat2*3 write(6,*) write(6,*)'Better Un FF file name:' read(5,'(a)')f3 write(6,*)'Number of atoms:' read(5,*)nat3 na3=3*nat3 write(6,*) write(6,*)'Better Per (resultant FF) file name:' read(5,'(a)')f4 nat4=nat1 na4=na1 C OPEN(20,FILE=f1,STATUS='OLD') CALL READFF(MX3,NA1,FCAR1) write(6,*)f1,' read in' CLOSE(20) OPEN(20,FILE=f2,STATUS='OLD') CALL READFF(MX3,NA2,FCAR2) write(6,*)f2,' read in' CLOSE(20) OPEN(20,FILE=f3,STATUS='OLD') CALL READFF(MX3,NA3,FCAR3) write(6,*)f3,' read in' CLOSE(20) c c diagonalized force field call TRED12(na1,FCAR1,b1,D1,2,IERR) write(6,*)'first FF diagonalized' c c diagonalize only part - na1 coordinates!: call TRED12(na1,FCAR2,b2,D2,2,IERR) write(6,*)'second FF diagonalized' c c diagonal(i,j)=b2(ii,i)*b2(jj,j)*FCAR(ii,jj)=d(i)*del(i,j) c c get scale factors write(6,*)'scale factors:' do 2 ii=1,na1 down=d1(ii) if (abs(down).lt.1.0d-5)then sf(ii)=1.0d0 else sf(ii)=d2(ii)/down endif write(6,6007)ii,sf(ii),d1(ii),d2(ii) 6007 format(i5,3f10.4) 2 continue c c diagonalize only part - na1 coordinates!: call TRED12(na1,FCAR3,b1,D1,2,IERR) write(6,*)'third FF diagonalized' c c scale diagonals: do 7 i=1,na1 7 D1(i)=D1(i)*sf(i) c c Back transformation: do 5 i=1,na1 do 5 j=1,na1 sum=0.0d0 do 6 ii=1,na1 6 sum=sum+b1(i,ii)*b1(j,ii)*D1(ii) 5 FCAR4(i,j)=sum c c rest of FF: do 8 i=na1+1,na3 do 8 j=1,na3 FCAR4(i,j)=FCAR3(i,j) 8 FCAR4(j,i)=FCAR3(j,i) ffac=0.0d0 gfac=0.0d0 do 1 i=1,na1 fac=sf(i) ffac=ffac+fac 1 gfac=gfac+fac**2 ffac=ffac/real(na1) gfac=gfac/real(na1) dev=sqrt(gfac-ffac**2) write(6,9000)ffac,dev 9000 format(' Average scale factor:',f7.4,' +/-',f7.4) c OPEN(20,FILE=f4) CALL WRITEFF(MX3,na1,fcar4) write(6,*)f4,' written' CLOSE(20) if(na3.gt.na1)then OPEN(20,FILE='COMPLETE.FC') CALL WRITEFF(MX3,na4,fcar4) write(6,*)f4,' COMPLETE.FC written' CLOSE(20) endif stop end C SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) C CONST=4.359828/0.5291772/0.5291772 CONST=1.0d0 N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) RETURN END SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) c CONST=4.359828/0.5291772**2 CONST=1.0d0 DO 6 I=1,N DO 6 J=1,N 6 FCAR(I,J)=FCAR(I,J)/CONST N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*300) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) COMMON/DIAG/E IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N DO 100 J = 1, I 100 Z(I,J) = A(I,J) 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) 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) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*300) DIMENSION Z(MX3,N),D(MX3),E(MX3) REAL*8 MACHEP COMMON/DIAG/E COMMON/INDEX/IND(MX3) 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 :::::::::: c 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 IY=IND(I) IND(I)=IND(K) IND(K)=IY 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