PROGRAM BUMATD c get diagonal intrinsic stretching constants IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,NAT0), 1B(MX,MX),NOB0,IS,UM(MX,MX) DIMENSION B0(MX,MX),LV(MX),IZ(NAT0) character*7 COO(MX),cty(7) character*2 as(MX) integer*4 d1(MX),d2(MX),d3(MX) data cty/'STRETCH','BEND ','XXXXXXX','TORSION','XXXXXXX', 1 'XXXXXXX','INTERMO'/ common/cl/d1,d2,d3,as,COO tol=1.0d-8 WRITE(*,*)' (SYMMETRIZED) B-MATRIX GENERATION' WRITE(*,*)' FILE.UMA ... INPUT FILE' OPEN(15,FILE='FILE.X') READ(15,*) READ(15,*)NOAT IF(NOAT.GT.NAT0)call report(' Too many atoms') do 21 I=1,NOAT 21 read(15,*)IZ(I) close(15) OPEN(15,FILE='FILE.UMA') READ(15,*) READ(15,*) READ(15,*)NOAT,NOB0 NA=NOAT*3 WRITE(*,*)NOAT,NOB0 ,MX if(NOB0.gt.MX)call report('Too many internal coordinates') DO 1001 I=1,NOB0 DO 1001 J=1,NA 1001 B(I,J)=0.0D0 if(NOB0.ne.NA-6)write(6,61000) 61000 format(/,' Warning: Number of Coordinates <> 3*NAT-6!',/) DO 20 I=1,NOAT READ(15,1020)(X(J,I),J=1,3),as(I) 1020 FORMAT(3F12.6,8x,a2) 20 CONTINUE IER=0 ISCAN=0 IS=0 NOB=0 277 NOB=NOB+1 if(NOB.gt.MX)call report('Too many internal coordinates') READ(15,*)ICODE,N1,N2,N3,N4,N5,N6,FACTOR write(6,1050)NOB,ICODE,N1,N2,N3,N4,N5,N6,FACTOR 1050 FORMAT(8I6,F12.6) if(ICODE.le.7)COO(NOB)=cty(ICODE) if(ICODE.eq.4)then d1(NOB)=n2 d2(NOB)=n3 d3(NOB)=0 else if(ICODE.eq.7)then if(NOB+5.gt.MX)call report('Too many internal coordinates') COO(NOB )='MSTRETC' COO(NOB+1)='MANGLE1' COO(NOB+2)='MANGLE2' COO(NOB+3)='MANGLE3' COO(NOB+4)='MANGLE4' COO(NOB+5)='MTORSIO' do i=NOB,NOB+5 d1(i)=n3 d2(i)=n4 enddo else if(ICODE.ge.8.and.ICODE.LE.13)then if(ICODE.eq. 8)COO(NOB)='MSTRETC' if(ICODE.eq. 9)COO(NOB)='ALPHA 1' if(ICODE.eq.10)COO(NOB)='ALPHA 2' if(ICODE.eq.11)COO(NOB)='BETA 1' if(ICODE.eq.12)COO(NOB)='BETA 2' if(ICODE.eq.13)COO(NOB)='TAU ' d1(NOB)=n3 d2(NOB)=n4 else d1(NOB)=n1 d2(NOB)=n2 d3(NOB)=n3 endif endif endif IF(DABS(FACTOR).LT.1.0D-6)FACTOR=1.0D0 IC=IABS(ICODE) IF(IC.EQ.0.OR.IC.GT.13) GO TO 220 GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13),IC 220 WRITE(*,1070)ICODE GO TO 200 1 CALL BOST GO TO 60 2 CALL BEND GO TO 60 3 CALL OPLA GO TO 60 4 CALL TORS GO TO 60 5 IS=IS+1 6 CALL LIBE GO TO 60 c intermolecular, together: 7 CALL INTER(IZ,ICODE) GO TO 60 c intermolecular, separate: 8 CALL INTER(IZ,ICODE) GO TO 60 9 CALL INTER(IZ,ICODE) GO TO 60 10 CALL INTER(IZ,ICODE) GO TO 60 11 CALL INTER(IZ,ICODE) GO TO 60 12 CALL INTER(IZ,ICODE) GO TO 60 13 CALL INTER(IZ,ICODE) 60 IF(IER.NE.0)GO TO 190 IF(DABS(FACTOR-1.0D0).LT.1.0D-8)GO TO 105 IF(IC.EQ.5)GO TO 80 70 ISW=0 I=NOB GO TO 90 80 ISW=1 I=NOB0+IS 90 DO 100 J=1,NA 100 B(I,J)=B(I,J)*FACTOR IF(ISW.EQ.1)GO TO 70 105 IF(ICODE.GT.0)GO TO 2775 C IF(ICODE.EQ.-5)GO TO 110 C ISW=0 C I=NOB-1 C GO TO 130 C 110 ISW=1 C 120 I=NOB-2 C 130 DO 140 J=1,NA C B(NOB0+IS,J) C 140 B(I,J)=B(I,J) + B(NOB,J) C NOB=NOB0+IS C ISW=0 C GO TO 120 2775 if(NOB.lt.NOB0)goto 277 write(6,*)'Coordinate definition read' DO 2770 I=1,NA DO 2770 J=1,NA 2770 UM(I,J)=0.0D0 READ(15,*)IUM IF(IUM.NE.0)THEN DO 2772 I=1,NA-6 READ(15,*)J SUM=0.0D0 DO 2774 IK=1,J READ(15,*)IJ,UM(I,IJ) LV(IK)=IJ 2774 SUM=SUM+UM(I,IJ)*UM(I,IJ) DO 2773 IK=1,J 2773 UM(I,LV(IK))=UM(I,LV(IK))/sqrt(SUM) 2772 CONTINUE WRITE(*,*)' USER-DEFINED U MATRIX' M=NA-6 ELSE DO 2771 I=1,NA 2771 UM(I,I)=1.0D0 M=NOB0 ENDIF CLOSE(15) IF(IUM.NE.0)THEN DO 1556 I=1,NA-6 DO 1556 J=1,NA SUM=0.0D0 DO 1557 K=1,NOB0 1557 SUM=SUM+UM(I,K)*B(K,J) 1556 B0(I,J)=SUM ELSE DO 1558 I=1,NOB0 DO 1558 J=1,NA 1558 B0(I,J)=B(I,J) ENDIF DO 155 I=1,NOB0 DO 155 J=1,NA 155 IF(DABS(B(I,J)).LT.tol) B(I,J) = 0.0D0 WRITE(*,*)' OUTPUT FILE: B.MAT' OPEN(9,FILE='B.MAT') WRITE(9,901)M 901 format(I6,' Internal coordinates') WRITE(9,903)NA 903 format(I6,' Cartesian coordinates') WRITE(9,902) 902 format(' Transformation matrix, non-zero elements:') do 15 I=1,M WRITE(9,9011)I 9011 format(I6,' internal:') nn=0 do 14 J=1,NA 14 if(DABS(B0(I,J)).gt.tol)nn=nn+1 write(9,9012)nn 9012 format(I6) nn=0 do 15 J=1,NA if(DABS(B0(I,J)).gt.tol)then nn=nn+1 ia=(J+2)/3 ix=J-3*(ia-1) WRITE(9,9100)ia,ix,B0(I,J) 9100 format(2I5,F15.8) endif 15 continue CLOSE(9) CALL CTI WRITE(*,*) WRITE(*,*)' >> PROGRAM TERMINATED <<' STOP 190 IF(IER.EQ.1)WRITE(*,1130) 200 IF(ISCAN.EQ.0)WRITE(16,1170) ISCAN=ISCAN+1 IER=0 1070 FORMAT(/,1X,'ILLEGAL CODE',I5,' CHOSEN') 1130 FORMAT(' ILLEGAL SPECIFICATION OF I, J, K, L, IX OR JX') 1170 FORMAT(/,1X,'*** PROGRAM WILL SCAN FOR FURTHER ERRORS IN DATA', $ ' DECK',/) END SUBROUTINE BOST C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A BOND STRETCH C AS DEFINED BY WILSON. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RIJ(3),RJK(3),RJL(3),EIJ(3),EJK(3),EKL(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 30 IF(J.LE.0.OR.J.GT.NOAT)GO TO 30 IF(K.NE.0)GO TO 30 IF(L.NE.0)GO TO 30 IF(IX.NE.0)GO TO 30 IF(JX.NE.0)GO TO 30 DIJSQ=0.D0 DO 10 M=1,3 T=X(M,J)-X(M,I) RIJ(M)=T 10 DIJSQ=DIJSQ+T*T DIJ=DSQRT(DIJSQ) II=3*(I-1) JJ=3*(J-1) DO 20 M=1,3 T=RIJ(M) IF(DABS(T).LT.1.D-8)GO TO 20 T=T/DIJ B(NOB,II+M)=-T B(NOB,JJ+M)=T 20 CONTINUE RETURN 30 IER=1 WRITE(*,200) 200 FORMAT('BOST ERROR') RETURN END SUBROUTINE BEND C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS OF A VALENCE C ANGLE BENDING COORDINATE AS DEFINED BY WILSON. C I AND K ARE THE NUMBERS OF THE END ATOMS. C J IS THE NUMBER OF THE CENTRAL ATOM. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EKL(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 50 IF(J.LE.0.OR.J.GT.NOAT)GO TO 50 IF(K.LE.0.OR.K.GT.NOAT)GO TO 50 IF(L.NE.0)GO TO 50 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 50 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 50 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DOTJ=0.D0 DO 30 M=1,3 T=RJI(M)/DJI EJI(M)=T TP=RJK(M)/DJK EJK(M)=TP 30 DOTJ=DOTJ+T*TP IF(DABS(DOTJ).GT.0.99995D0)GO TO 60 SINJ=DSQRT(1.D0-DOTJ*DOTJ) II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) DO 40 M=1,3 SMI=DX*(DOTJ*EJI(M)-EJK(M))/(DJI*SINJ) IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=SMI SMK=DX*(DOTJ*EJK(M)-EJI(M))/(DJK*SINJ) IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=SMK SUM=SMI+SMK 40 IF(DABS(SUM).GE.1.D-8)B(NOB,JJ+M)=-SUM RETURN 50 IER=1 *C*C* ADDED BY KATH SPRING 1991 TO FIND INPUT ERROR * NEXT TWO LINES WRITE(*,200) 200 FORMAT('BEND ERROR') RETURN 60 IER=-1 WRITE(*,1000) RETURN C 1000 FORMAT(' I-J-K IS COLINEAR - USE LINEAR BEND') C END SUBROUTINE OPLA C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR AN OUT OF C PLANE WAGGING COORDINATE AS DEFINED BY DECIUS, MCINTOSH, MICHAELIAN C AND PETERSON. SUBROUTINE CODED BY M PETERSON, UNIV OF TORONTO. C I IS THE END ATOM (ATOM WAGGED WITH RESPECT TO J-K-L PLANE). C J IS THE APEX ATOM (ATOMS I, K AND L ARE ATTACHED TO J). C K AND L ARE THE ANCHOR ATOMS (DEFINE THE J-K-L PLANE). C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EJL(3) DIMENSION C1(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 60 IF(J.LE.0.OR.J.GT.NOAT)GO TO 60 IF(K.LE.0.OR.K.GT.NOAT)GO TO 60 IF(L.LE.0.OR.L.GT.NOAT)GO TO 60 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 60 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 60 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DJLSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,L)-X(M,J) RJL(M)=T DJLSQ=DJLSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DJL=DSQRT(DJLSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 COSI=0.D0 COSK=0.D0 COSL=0.D0 DO 30 M=1,3 T=RJI(M)/DJI EJI(M)=T TP=RJK(M)/DJK EJK(M)=TP TPP=RJL(M)/DJL EJL(M)=TPP COSI=COSI+TP*TPP COSK=COSK+T*TPP 30 COSL=COSL+T*TP IF(DABS(COSI).GT.0.99995D0)GO TO 70 SINSIN=1.D0-COSI*COSI SINI=DSQRT(SINSIN) C1(1)=EJK(2)*EJL(3)-EJK(3)*EJL(2) C1(2)=EJK(3)*EJL(1)-EJK(1)*EJL(3) C1(3)=EJK(1)*EJL(2)-EJK(2)*EJL(1) DOT=EJI(1)*C1(1)+EJI(2)*C1(2)+EJI(3)*C1(3) SINT=DOT/SINI IF(DABS(SINT).GT.0.00005D0)WRITE(*,1020) IF(DABS(SINT).GT.0.99995D0)GO TO 80 COST=DSQRT(1.D0-SINT*SINT) TANT=SINT/COST II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) LL=3*(L-1) COSSIN=COST*SINI DO 50 M=1,3 T=C1(M)/COSSIN SMI=(T-TANT*EJI(M))/DJI IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=DX*SMI SMK=T*(COSI*COSK-COSL)/(SINSIN*DJK) IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=DX*SMK SML=T*(COSI*COSL-COSK)/(SINSIN*DJL) IF(DABS(SML).GE.1.D-8)B(NOB,LL+M)=DX*SML SUM=SMI+SMK+SML 50 IF(DABS(SUM).GE.1.D-8)B(NOB,JJ+M)=-DX*SUM RETURN 60 IER=1 *C*C ADDED BY KATH, NEXT TWO LINES, SPRING 1991, LOOKING FOR ERRORS * WRITE(*,200) 200 FORMAT('OPLA ERROR') RETURN 70 IER=-1 WRITE(*,1000) RETURN 80 IER=-1 WRITE(*,1010) RETURN C 1000 FORMAT(' K-J-L IS COLINEAR (NO PLANE DEFINED FOR WAG OF I)') 1010 FORMAT(' I IS PERPENDICULAR TO J-K-L PLANE - USE VALENCE ANGLE 1 BENDS') 1020 FORMAT('+',86X,'*** WARNING: WAG OF A NON-PLANAR SYSTEM ***') C END SUBROUTINE TORS C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR TORSION AS DEFIN C BY R L HILDERBRANDT IN J MOLEC SPEC, 44, 599 (1972). C SUBROUTINE CODED BY M PETERSON, DEPT OF CHEMISTRY, UNIV OF TORONTO. C J AND K DEFINE THE BOND UNDER TORSION. C NI AND NL ARE THE NUMBER OF ATOMS ATTACHED TO J AND K RESPECTIVELY C (NI, NL <= 5). 2 DATA CARDS ARE READ: (1) CONTAINS NI ATOM NUMBERS C FOR THE I-TYPE ATOMS, AND (2) CONTAINS NL ATOM NUMBERS FOR THE L-TY C ATOMS (BOTH CARDS ARE IN 5I4 FORMAT). C C IATOM, LATOM: ATOM NUMBERS FOR THE I- AND L-TYPE ATOMS (SIZE: 5) C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,NI,J,K,NL,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RIJ(3),RJK(3),RLK(3),EIJ(3),EJK(3),ELK(3) DIMENSION CR(3),IATOM(5),LATOM(5),SJ(3),SK(3) C READ(15,*)(IATOM(I),I=1,NI) READ(15,*)(LATOM(L),L=1,NL) IF(NI.LE.0.OR.NI.GT.5)GO TO 110 IF(J.LE.0.OR.J.GT.NOAT)GO TO 110 IF(K.LE.0.OR.K.GT.NOAT)GO TO 110 IF(NL.LE.0.OR.NL.GT.5)GO TO 110 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 110 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 110 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJKSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 SJ(M)=0.D0 SK(M)=0.D0 T=X(M,K)-X(M,J) RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJK=1.D0/DSQRT(DJKSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DO 30 M=1,3 30 EJK(M)=RJK(M)*DJK JJ=3*(J-1) KK=3*(K-1) C C LOOP OVER THE I-TYPE ATOMS C DO 60 N=1,NI I=IATOM(N) IF(I.LE.0.OR.I.GT.NOAT)GO TO 110 DIJSQ=0.D0 DO 40 M=1,3 T=X(M,J)-X(M,I) RIJ(M)=T 40 DIJSQ=DIJSQ+T*T DIJ=1.D0/DSQRT(DIJSQ) COSJ=0.D0 DO 50 M=1,3 T=RIJ(M)*DIJ EIJ(M)=T 50 COSJ=COSJ-T*EJK(M) IF(DABS(COSJ).GT.0.99995D0)GO TO 120 SIN2J=(1.D0-COSJ*COSJ)*DFLOAT(NI) II=3*(I-1) CR(1)=EIJ(2)*EJK(3)-EIJ(3)*EJK(2) CR(2)=EIJ(3)*EJK(1)-EIJ(1)*EJK(3) CR(3)=EIJ(1)*EJK(2)-EIJ(2)*EJK(1) DO 60 M=1,3 T=CR(M)/SIN2J SMI=T*DIJ IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=-DX*SMI SMK=T*COSJ*DJK SK(M)=SK(M)+SMK SMJ=SMI-SMK 60 SJ(M)=SJ(M)+SMJ C C LOOP OVER THE L-TYPE ATOMS C DO 90 N=1,NL L=LATOM(N) IF(L.LE.0.OR.L.GT.NOAT)GO TO 110 DLKSQ=0.D0 DO 70 M=1,3 T=X(M,K)-X(M,L) RLK(M)=T 70 DLKSQ=DLKSQ+T*T DLK=1.D0/DSQRT(DLKSQ) COSK=0.D0 DO 80 M=1,3 T=RLK(M)*DLK ELK(M)=T 80 COSK=COSK+EJK(M)*T IF(DABS(COSK).GT.0.99995D0)GO TO 120 SIN2K=(1.D0-COSK*COSK)*DFLOAT(NL) LL=3*(L-1) CR(1)=ELK(3)*EJK(2)-ELK(2)*EJK(3) CR(2)=ELK(1)*EJK(3)-ELK(3)*EJK(1) CR(3)=ELK(2)*EJK(1)-ELK(1)*EJK(2) DO 90 M=1,3 T=CR(M)/SIN2K SML=T*DLK IF(DABS(SML).GE.1.D-8)B(NOB,LL+M)=-DX*SML SMJ=T*COSK*DJK SJ(M)=SJ(M)+SMJ SMK=SML-SMJ 90 SK(M)=SK(M)+SMK DO 100 M=1,3 SMJ=SJ(M) IF(DABS(SMJ).GE.1.D-8)B(NOB,JJ+M)=SMJ*DX SMK=SK(M) 100 IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=SMK*DX RETURN 110 IER=1 WRITE(*,200) I,L 200 FORMAT('TORS ERROR AT I=',I2,' AND L=',I2) RETURN 120 IER=-1 WRITE(*,1030) RETURN 1030 FORMAT(' I-J-K OR J-K-L IS COLINEAR (NO TORSION POSSIBLE)') END SUBROUTINE LIBE C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A LINEAR BEND C OR FOR A PAIR OF PERPENDICULAR LINEAR BENDS. C I AND K ARE THE END ATOMS. C J IS THE CENTRAL ATOM. C C A GIVES THE CARTESIAN COORDINATES OF A POINT IN SPACE, SUCH C THAT THE VECTOR FROM ATOM J TO POINT A IS PERPENDICULAR TO C THE LINE I-J-K AND SERVES TO ORIENT THE COORDINATES IN SPACE. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),EJK(3),UP(3),UN(3),UNIT(3) DIMENSION A(3) C READ(15,1000)A IF(I.LE.0.OR.I.GT.NOAT)GO TO 60 IF(J.LE.0.OR.J.GT.NOAT)GO TO 60 IF(K.LE.0.OR.K.GT.NOAT)GO TO 60 IF(L.NE.0)GO TO 60 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 60 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 60 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DXSQ=0.D0 DJASQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) DXSQ=DXSQ+T*T T=A(M)-TP UN(M)=T 20 DJASQ=DJASQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DX=DSQRT(DXSQ) DJA=DSQRT(DJASQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DOTJ=0.D0 DOTP=0.D0 DO 30 M=1,3 T=RJI(M)/DJI TP=RJK(M)/DJK EJK(M)=TP DOTJ=DOTJ+T*TP TP=UN(M)/DJA UNIT(M)=TP 30 DOTP=DOTP+T*TP TEST=DABS(DOTJ)-1.D0 IF(DABS(TEST).GT.0.00005D0)GO TO 70 IF(DABS(DOTP).GT.0.00005D0)GO TO 80 II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) DO 40 M=1,3 T=UNIT(M) IF(DABS(T).LT.1.D-8)GO TO 40 T=-DX*T SMI=T/DJI B(NOB,II+M)=SMI SMK=T/DJK B(NOB,KK+M)=SMK B(NOB,JJ+M)=-SMI-SMK 40 CONTINUE IF(IC.EQ.6)RETURN NOB=NOB0+IS UP(1)=EJK(2)*UNIT(3)-EJK(3)*UNIT(2) UP(2)=EJK(3)*UNIT(1)-EJK(1)*UNIT(3) UP(3)=EJK(1)*UNIT(2)-EJK(2)*UNIT(1) DO 50 M=1,3 T=UP(M) IF(DABS(T).LT.1.D-8)GO TO 50 T=-DX*T SMI=T/DJI B(NOB,II+M)=SMI SMK=T/DJK B(NOB,KK+M)=SMK B(NOB,JJ+M)=-SMI-SMK 50 CONTINUE RETURN 60 IER=1 RETURN 70 IER=-1 WRITE(*,1020) RETURN 80 IER=-1 WRITE(*,1030) RETURN 1000 FORMAT(3D23.16) 1020 FORMAT(' I-J-K NOT COLINEAR - USE VALENCE ANGLE BEND') 1030 FORMAT(' ATOM A NOT PERPENDICULAR TO I-J-K AT J') END SUBROUTINE CTI IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0,NX=MX) DIMENSION FINT(MX,MX),B(MX,NX),FCAR(MX,MX),nb(MX),lb(MX,NX) allocatable::am(:,:),ami(:,:),cm(:),t(:,:) character*7 COO(MX) character*2 as(MX) integer*4 d1(MX),d2(MX),d3(MX) common/cl/d1,d2,d3,as,COO WRITE(*,*)' INTRISIC FORCE FIELD GENERATION' open(9,file='B.MAT') read(9,*)nob0 read(9,*)na read(9,*) do 13 i=1,nob0 read(9,*) read(9,*)nb(i) if(nb(i).gt.NX)call report('too many B-matrix elements') do 13 ii=1,nb(i) read(9,*)ia,ix,b(i,ii) 13 lb(i,ii)=ix+3*(ia-1) close(9) write(6,*)'B.MAT read, ',na,' atomic ',nob0,' internal' NAT3=NA WRITE(*,*)' FILE.FC ... INPUT FILE WITH' WRITE(*,*)' THE CARTESIAN FORCE FIELD' OPEN(20,FILE='FILE.FC',STATUS='OLD') N=NAT3 NTT= -4 NTE= 5 IF (NTE.GT.N) GO TO 120 110 NTT=NTT+5 DO 100 I=NTT,N NT=I IF (NT.GT.NTE) NT=NTE 100 READ(20,27) (FCAR(I,J),J=NTT,NT) 27 FORMAT(4X,5D14.6) NTE=NTE+5 IF (NTE.LE.N) GOTO 110 120 NTT=NTT+5 DO 130 I=NTT,N 130 READ(20,27) (FCAR(I,J),J=NTT,I) CLOSE(20) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) DO 6 I=1,N DO 6 J=1,N C 6 FCAR(I,J)=FCAR(I,J)* 8.23D0 / 0.527D0 6 FCAR(I,J)=FCAR(I,J)*4.3597482D0/(0.529177249D0*0.529177249D0) c 1Hartree=4.3597482e-18 J; 1Bohr=0.529177249A write(6,*)'FILE.FC read' TOL= 1.0D-33 allocate(am(nob0,nob0),ami(nob0,nob0),t(nob0,nob0)) write(6,*)'fields allocated' alpha=1.0d-6 c c do 51 i=1,nob0 do 51 j=1,nob0 sum=0.0d0 do 511 lc=1,nb(i) do 511 kc=1,nb(j) 511 if(lb(i,lc).eq.lb(j,kc))sum=sum+b(i,lc)*b(j,kc) 51 t(i,j)=sum write(6,*)'t made' do 5 i=1,nob0 do 5 j=1,nob0 5 am(i,j)=t(i,j)**2 do 52 i=1,nob0 52 am(i,i)=am(i,i)+alpha write(6,*)'a done' CALL INV(am,ami,nob0,TOL) deallocate(am) allocate(cm(nob0)) do 7 i=1,nob0 cm(i)=0.0d0 do 7 lc=1,nb(i) do 7 ic=1,nb(i) 7 cm(i)=cm(i)+b(i,ic)*b(i,lc)*FCAR(lb(i,ic),lb(i,lc)) write(6,*)'c done' do 4 i=1,nob0 do 4 j=1,nob0 FINT(i,j)=0 if(i.eq.j)then do 41 ip=1,nob0 41 FINT(i,j)=FINT(i,j)+ami(i,ip)*cm(ip) endif 4 continue deallocate(cm,t) deallocate(ami) OPEN(2,FILE='INTY.FC') WRITE(2,*)' Force constants in internal coordinates:' CALL WRITEFF(MX,nob0,FINT) close(2) WRITE(*,*)'INTY.FC : OUTPUT WITH THE INTRISIC FORCE CONSTANTS' OPEN(2,FILE='IFC.TXT') do 8 i=1,nob0 8 write(2,800)i,COO(i),d1(i),d2(i),FINT(i,i) 800 format(i6,' ',a7,i6,i6,f12.6) close(2) RETURN END SUBROUTINE INV(MA,IA,N,TOL) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) REAL*8 MA,IA INTEGER*4 oo,n2 DIMENSION IA(N,N),MA(N,N) real*8, allocatable ::e(:) write(6,*)N*2*N/8,' bytes',N,2*N allocate(e(N*2*N)) n2=2*N write(6,*)'Inversion started' DO 1 ii=1,N DO 1 jj=1,N e(n2*(ii-1)+jj)=ma(ii,jj) e(n2*(ii-1)+jj+n)=0.0D0 1 if (ii.EQ.jj)e(n2*(ii-1)+jj+n)=1.0D0 write(6,*)'E initialized' DO 2 ii=1,n-1 IW=II if (DABS(e(n2*(iw-1)+iw)).LT.TOL) then DO 3 io=iw+1,n oo=IO if (io.GT.n)oo=io-n if (DABS(e(n2*(oo-1)+iw)).GE.TOL) goto 11 3 CONTINUE do 13 i=1,N 13 write(6,90)(MA(i,j),j=1,N) 90 format(12g12.4) WRITE(*,*)' INVERSE MATRIX CAN NOT BE FOUND ' STOP 11 CONTINUE DO 4 kk=1, 2*n w=e(n2*(iw-1)+kk) e(n2*(iw-1)+kk)=e(n2*(oo-1)+kk) 4 e(n2*(oo-1)+kk)=w ENDIF DO 5 jj=ii+1,n af=e(n2*(jj-1)+ii)/e(n2*(ii-1)+ii) DO 6 kk=ii+1, 2*n 6 e(n2*(jj-1)+kk)=e(n2*(jj-1)+kk)-e(n2*(ii-1)+kk)*af 5 e(n2*(jj-1)+ii)=0.0D0 2 CONTINUE write(6,*)'Inversion first loop' DO 7 ii=1,n-1 i2=n-ii+1 DO 8 jj=1,i2-1 j2=i2-jj af=e(n2*(j2-1)+i2)/e(n2*(i2-1)+i2) DO 9 kk=1, n 9 e(n2*(j2-1)+kk+n)=e(n2*(j2-1)+kk+n)-e(n2*(i2-1)+kk+n)*af e(n2*(j2-1)+i2)=0.0D0 8 CONTINUE 7 CONTINUE write(6,*)'Inversion second loop' DO 10 ii=1,n DO 12 jj=1,n 12 iA(ii,jj)=e(n2*(ii-1)+jj+n)/e(n2*(ii-1)+ii) 10 CONTINUE deallocate(E) RETURN END SUBROUTINE WRITEFF(MX,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX,MX) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N write(2,18)(ij,ij=N1,N3) 18 FORMAT(4x,5I14) DO 130 LN=N1,N 130 WRITE(2,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 report(s) character*(*) s write(6,*)s stop end SUBROUTINE INTER(IZ,ICODE) c coordinates defining 6 inter-molecular coordinates c c (list1) (list2) c atoms 1....N1 atoms 1..N2 c (molecule 1) (molecule 2) c implicit none integer*4 MX,IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,NOB0,IS,I,ICODE, 1IZ(*),MENDELEV,ix,jx,ia,IDIF,IERR,II,is1,is2,KA,KK,KX,NAT0 PARAMETER (NAT0=1000,MX=3*NAT0,MENDELEV=89) real*8 X,B,UM,th1(3,3),th2(3,3),c1(3),c2(3),m1,m2,xi,xj,c0(3), 1u1(3,3),u2(3,3),d1(3),d2(3),PI,STEP,FF,acc, 1v1(3),v2(3),e1(3),e2(3),D,T,A1,A2,B1,B2,D0,T0,A10,A20,B10,B20 COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) integer*4,allocatable:: list1(:),list2(:) real*4 AM(MENDELEV) data AM/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ PI=4.0d0*atan(1.0d0) allocate(list1(N1),list2(N2)) C READ(15,*)(list1(I),I=1,N1) READ(15,*)(list2(I),I=1,N2) m1=0.0d0 do 12 i=1,N1 12 m1=m1+dble(AM(IZ(list1(i)))) m2=0.0d0 do 13 i=1,N2 13 m2=m2+dble(AM(IZ(list2(i)))) c c numerical derivatives STEP=0.01D0 do 100 IDIF=0,6*NOAT if(IDIF.gt.0)THEN if(IDIF.GT.3*NOAT)THEN KK=IDIF-3*NOAT FF=-1.0D0 ELSE KK=IDIF FF=1.0D0 ENDIF KA=(KK+2)/3 KX=KK-3*(KA-1) X(KX,KA)=X(KX,KA)+FF*STEP ELSE FF=0.0d0 KA=0 KX=0 ENDIF c calculate centers of mass: do 3 ix=1,3 c1(ix)=0.0d0 do 4 i=1,N1 4 c1(ix)=c1(ix)+X(ix,list1(i))*dble(AM(IZ(list1(i)))) c1(ix)=c1(ix)/m1 c2(ix)=0.0d0 do 5 i=1,N2 5 c2(ix)=c2(ix)+X(ix,list2(i))*dble(AM(IZ(list2(i)))) 3 c2(ix)=c2(ix)/m2 d=0.0d0 do 6 ix=1,3 c0(ix)=c1(ix)-c2(ix) 6 d=d+c0(ix)**2 d=dsqrt(d) do 7 ix=1,3 7 c0(ix)=c0(ix)/d c c inter-molecular stretching - analytically: if(IDIF.EQ.0)then do 8 ia=1,NOAT c is atom ia part of the first or second molecule?: is1=0 is2=0 do 9 i=1,N1 9 if(list1(i).eq.ia)is1=1 do 10 i=1,N2 10 if(list2(i).eq.ia)is2=1 if(is1.ne.0.or.is2.ne.0)then do 11 ix=1,3 ii=ix+3*(ia-1) if(ICODE.EQ.7.or.ICODE.EQ.8)then if(is1.eq.1)B(NOB,ii)= c0(ix)*dble(AM(IZ(ia)))/m1 if(is2.eq.1)B(NOB,ii)=-c0(ix)*dble(AM(IZ(ia)))/m2 endif 11 continue endif 8 continue endif c calculate moments of inertia: do 1 ix=1,3 do 1 jx=1,3 th1(ix,jx)=0.0d0 th2(ix,jx)=0.0d0 do 2 i=1,N1 ia=list1(i) xi=x(ix,ia)-c1(ix) xj=x(jx,ia)-c1(jx) 2 th1(ix,jx)=th1(ix,jx)+xi*xj*dble(AM(IZ(ia))) do 1 i=1,N2 ia=list2(i) xi=x(ix,ia)-c2(ix) xj=x(jx,ia)-c2(jx) 1 th2(ix,jx)=th2(ix,jx)+xi*xj*dble(AM(IZ(ia))) c diagonalize moments of inertia: CALL TRED12(3,th1,u1,d1,2,IERR) CALL TRED12(3,th2,u2,d2,2,IERR) if(IDIF.EQ.0)write(6,607)d1,d2 607 format(' Inertias:',6F10.3) c angles between longest moment axes and the distance vecto A1=acc(c0(1)*u1(1,1)+c0(2)*u1(2,1)+c0(3)*u1(3,1)) A2=acc(c0(1)*u2(1,1)+c0(2)*u2(2,1)+c0(3)*u2(3,1)) c c arbitrary vectors defining rotating axes v1(1)=u1(2,1)*c0(3)-u1(3,1)*c0(2) v1(2)=u1(3,1)*c0(1)-u1(1,1)*c0(3) v1(3)=u1(1,1)*c0(2)-u1(2,1)*c0(1) v2(1)=u2(2,1)*c0(3)-u2(3,1)*c0(2) v2(2)=u2(3,1)*c0(1)-u2(1,1)*c0(3) v2(3)=u2(1,1)*c0(2)-u2(2,1)*c0(1) e1(1)=v1(2)*u1(3,1)-v1(3)*u1(2,1) e1(2)=v1(3)*u1(1,1)-v1(1)*u1(3,1) e1(3)=v1(1)*u1(2,1)-v1(2)*u1(1,1) e2(1)=v2(2)*u2(3,1)-v2(3)*u2(2,1) e2(2)=v2(3)*u2(1,1)-v2(1)*u2(3,1) e2(3)=v2(1)*u2(2,1)-v2(2)*u2(1,1) c angles between the second longest u and the C0U plane B1=acc(e1(1)*u1(1,2)+e1(2)*u1(2,2)+e1(3)*u1(3,2)) B2=acc(e2(1)*u2(1,2)+e2(2)*u2(2,2)+e2(3)*u2(3,2)) c torsion angle T=acc(v1(1)*v2(2)+v1(2)*v2(2)+v1(3)*v2(3)) if(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3).lt.0.0d0)T=-T c try to keep conventions: c torsion -180...180 if(T.lt.-180.0d0)T= 360.0d0+T if(T.gt.+180.0d0)T=-360.0d0+T c axis angles: the sharp ones if(A1.gt.90.0d0)A1=180.0d0-A1 if(A2.gt.90.0d0)A2=180.0d0-A2 if(B1.gt.90.0d0)B1=180.0d0-B1 if(B2.gt.90.0d0)B2=180.0d0-B2 c write(6,609)D,A1,A2,B1,B2,T if(IDIF.eq.0)then D0=d A10=A1 A20=A2 B10=B1 B20=B2 T0=1 else c B(NOB ,KK)=B(NOB ,KK)+FF*d if(ICODE.EQ.7)then B(NOB+1,KK)=B(NOB+1,KK)+FF*A1 B(NOB+2,KK)=B(NOB+2,KK)+FF*A2 B(NOB+3,KK)=B(NOB+3,KK)+FF*B1 B(NOB+4,KK)=B(NOB+4,KK)+FF*B2 B(NOB+5,KK)=B(NOB+5,KK)+FF*T else if(ICODE.EQ. 9)B(NOB,KK)=B(NOB,KK)+FF*A1 if(ICODE.EQ.10)B(NOB,KK)=B(NOB,KK)+FF*A2 if(ICODE.EQ.11)B(NOB,KK)=B(NOB,KK)+FF*B1 if(ICODE.EQ.12)B(NOB,KK)=B(NOB,KK)+FF*B2 if(ICODE.EQ.13)B(NOB,KK)=B(NOB,KK)+FF*T endif ENDIF if(IDIF.NE.0)X(KX,KA)=X(KX,KA)-FF*STEP 100 continue if(ICODE.EQ.7)then do 101 i=1,5 do 101 KK=1,3*NOAT 101 B(NOB+i,KK)=B(NOB+i,KK)/2.0d0/STEP*PI/180.0d0 NOB=NOB+5 write(6,609)D0,A10,A20,B10,B20,T0,NOB-6,NOB 609 format(' D A1 A2 B1 B2 T:',F10.4,5F8.2,i4,' -',i4) else if(ICODE.GT.8)then do 102 KK=1,3*NOAT 102 B(NOB,KK)=B(NOB,KK)/2.0d0/STEP*PI/180.0d0 endif write(6,609)D0,A10,A20,B10,B20,T0,NOB endif 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 function acc(p) implicit none real*8 A,acc,p,PI PI=4.0d0*atan(1.0d0) A=0.0D0 IF(p.LE.-1.0D0)A=PI IF(DABS(p).LT.1.0D0)A=ACOS(p) acc=180.0D0*A/PI return end