PROGRAM BMATRIX IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (M0=3*145) CHARACTER*20 TITLE logical lex COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,M0/3), 1B(M0,M0),NOB0,IS,UM(M0,M0) DIMENSION Q(M0/3), 2Z(M0/3),CM(3),R(M0/3,3),BI(M0,M0),B0(M0,M0) WRITE(6,*)' B-MATRIX GENERATION' WRITE(6,*)' FILE.UMA ... INPUT FILE' DO 10011 I=1,M0 DO 10011 J=1,M0 10011 B(I,J)=0.0D0 OPEN(15,FILE='FILE.UMA') READ(15,1001)TITLE READ(15,1001)TITLE 1001 FORMAT(A20) READ(15,*)NOAT,NOB0 WRITE(6,*)TITLE NA=NOAT*3 DO 20 I=1,NOAT READ(15,1020)(X(J,I),J=1,3) 1020 FORMAT(3F12.6) 20 CONTINUE IER=0 ISCAN=0 IS=0 DO 277 NOB=1,NOB0 READ(15,1050)ICODE,N1,N2,N3,N4,N5,N6,FACTOR 1050 FORMAT(7I4,F12.6) IF(DABS(FACTOR).LT.1.0D-6)FACTOR=1.0D0 IC=IABS(ICODE) IF(IC.EQ.0.OR.IC.GT.6) GO TO 220 GO TO (1,2,3,4,5,6),IC 220 WRITE(6,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 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 277 277 CONTINUE DO 2770 I=1,M0 DO 2770 J=1,M0 2770 UM(I,J)=0.0D0 READ(15,*)IUM IF(IUM.NE.0)THEN WRITE(6,*)' USER-DEFINED U MATRIX' DO 2772 I=1,NA-6 READ(15,*)J DO 2772 IK=1,J 2772 READ(15,*)IJ,UM(I,IJ) ELSE DO 2771 I=1,M0 2771 UM(I,I)=1.0D0 ENDIF CLOSE(15) DO 155 I=1,NOB0 DO 155 J=1,NA IF(DABS(B(I,J)).LT.1.0D-08) B(I,J) = 0.0D0 155 CONTINUE 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,NA-6 DO 1558 J=1,NA 1558 B0(I,J)=B(I,J) ENDIF WRITE(6,*)' Redundant coordinate generation ...' OPEN(15,FILE='FTRY.INP') READ(15,*) READ(15,*) NOSC DO 1041 I=1,NOSC 1041 READ(15,*) READ(15,*)(Q(I),I=1,NOAT) 1040 FORMAT(6G12.6) CHATOT=0.0 NOAT=NA/3 DO 77 I=1,NOAT 77 CHATOT=CHATOT+Q(I) WRITE(6,777)CHATOT 777 FORMAT(' Total CHARGE',F10.6) 778 FORMAT(' Total MASS ',F10.6) C Read-in atomic masses: READ(15,*)NM IF(NM.NE.1)THEN WRITE(6,*)' Only one molecule allowed here !' CLOSE(15) STOP ENDIF READ(15,1040)(Z(I),I=1,NOAT) CLOSE(15) BM=0.0 DO 102 IA=1,NOAT 102 BM=BM+Z(IA) WRITE(6,778)BM C Read-in coordinates DO 107 I=1,3 107 CM(I)=0.0 inquire(FILE='FILE.X',exist=lex) if(lex)then OPEN(15,FILE='FILE.X') READ(15,*) READ(15,*)natx write(6,*)natx,' atoms in FILE.X' else OPEN(15,FILE='FILE.XYZ',status='old') READ(15,*) endif DO 106 IA=1,NOAT READ(15,*)IDUM,(R(IA,IX),IX=1,3) DO 106 IX=1,3 106 CM(IX)=CM(IX)+Z(IA)*R(IA,IX)/BM DO 108 IA=1,NOAT DO 108 IX=1,3 108 R(IA,IX)=R(IA,IX)-CM(IX) CLOSE(15) c C Define redundant coordinates: BMP=1.0/SQRT(BM) DO 101 IA=1,3 C Translation along Xia ICO=NOB0+IA DO 103 IX=1,3*NOAT 103 B0(ICO,IX)=0.0 DO 101 I=1,NOAT 101 B0(ICO,3*I-3+IA)=BMP*Z(I) DO 104 IA=1,3 C Rotation about Xia ICO=NOB0+3+IA DO 110 IX=1,3*NOAT 110 B0(ICO,IX)=0.0 IB=IA+1 IF(IB.GT.3)IB=1 IC=IB+1 IF(IC.GT.3)IC=1 AII=0.0 DO 205 I=1,NOAT 205 AII=AII+Z(I)*(R(I,IB)**2+R(I,IC)**2) AIIB=1.0/SQRT(AII) DO 104 I=1,NOAT B0(ICO,3*I-3+IB)=-AIIB*R(I,IC)*Z(I) 104 B0(ICO,3*I-3+IC)= AIIB*R(I,IB)*Z(I) C WRITE(6,*)' B DONE' WRITE(6,*)' OUTPUT FILE: B.MAT' TOL=0.00000001 CALL INV(B0,BI,NA,TOL) WRITE(6,*)' Inversed B-matrix calculated' OPEN(9,FILE='B.MAT',FORM='UNFORMATTED') WRITE(9)NOB0,NA,((B0(I,J),J=1,NA),I=1,NA), 1 ((BI(I,J),J=1,NA),I=1,NA) CLOSE(9) WRITE(6,*)' >> PROGRAM TERMINATED <<' STOP 190 IF(IER.EQ.1)WRITE(6,1130) 200 IF(ISCAN.EQ.0)WRITE(16,1170) ISCAN=ISCAN+1 IER=0 WRITE(6,1120) 1070 FORMAT(/,1X,'ILLEGAL CODE',I5,' CHOSEN') 1120 FORMAT(////,1X,'*** END OF B MATRIX COMPUTATION',//) 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 (M0=3*145) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,M0/3),B(M0,M0) 1,NOB0,IS,UM(M0,M0) 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(6,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 (M0=3*145) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,M0/3),B(M0,M0) 1,NOB0,IS,UM(M0,M0) 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(6,200) 200 FORMAT('BEND ERROR') RETURN 60 IER=-1 WRITE(6,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 (M0=3*145) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,M0/3),B(M0,M0) 1,NOB0,IS,UM(M0,M0) 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(6,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(6,200) 200 FORMAT('OPLA ERROR') RETURN 70 IER=-1 WRITE(6,1000) RETURN 80 IER=-1 WRITE(6,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 (M0=3*145) COMMON/BMAT/IC,NI,J,K,NL,IX,JX,NOAT,NOB,IER,X(3,M0/3),B(M0,M0) 1,NOB0,IS,UM(M0,M0) 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,1000)(IATOM(I),I=1,NI) READ(15,1000)(LATOM(L),L=1,NL) 1000 FORMAT(16I4) 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(6,200) I,L 200 FORMAT('TORS ERROR AT I=',I2,' AND L=',I2) RETURN 120 IER=-1 WRITE(6,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 (M0=3*145) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,M0/3),B(M0,M0) 1,NOB0,IS,UM(M0,M0) 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(6,1020) RETURN 80 IER=-1 WRITE(6,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 INV(MA,IA,N,TOL) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (M0=3*145) REAL*8 MA,IA INTEGER OO DIMENSION IA(M0,M0),MA(M0,M0),E(M0,M0*2) NMAT=N DO 1 ii=1,nmat DO 1 jj=1,nmat e(ii,jj)=ma(ii,jj) E(II,JJ+NMAT)=0.0D0 if (ii.EQ.jj)e(ii,jj+nmat)=1.0D0 1 CONTINUE DO 2 ii=1,nmat-1 IW=II if (DABS(e(iw,iw)).LT.TOL) then DO 3 io=iw+1,nmat OO=IO if (io.GT.nmat) oo=io-nmat if (DABS(e(oo,iw)).GE.TOL) goto 11 3 CONTINUE WRITE(6,*)' INVERSE MATRIX CAN NOT BE FOUND ' STOP 11 CONTINUE DO 4 kk=1, 2*nmat w=e(iw,kk) e(iw,kk)=e(oo,kk) 4 e(oo,kk)=w ENDIF DO 5 jj=ii+1,nmat DO 6 kk=ii+1, 2*nmat 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii) 5 e(jj,ii)=0.0D0 2 CONTINUE DO 7 ii=1, nmat-1 i2=nmat-ii+1 DO 8 jj=1,i2-1 j2=i2-jj DO 9 kk=1, nmat 9 e(j2,kk+nmat)=e(j2,kk+nmat)- 1e(i2,kk+nmat)*e(j2,i2)/e(i2,i2) e(j2,i2)=0 8 CONTINUE 7 CONTINUE DO 10 ii=1,nmat DO 12 jj=1,nmat 12 iA(ii,jj)=e(ii,jj+nmat)/e(ii,ii) 10 CONTINUE RETURN END